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I.  INTRODUCTION 


In  recent  years,  two  developing  concepts  have  been  considered  by  mission 
designers  to  reduce  the  required  propellant  mass  for  interplanetary  missions,  thus 
allowing  for  increased  mass  of  scientific  payloads.  Low  thrust  propulsion  allows  for 
greatly  improved  fuel  efficiency  and  has  the  potential  to  increase  payload  mass  fractions 
as  well  as  providing  trajectories  not  possible  with  impulsive  thrust.  Aerocapture  is  the 
careful  management  of  a  hypersonic  atmospheric  pass  prior  to  orbit  insertion  to  remove 
excess  arrival  velocity  of  a  spacecraft  obviating  the  need  for  a  large  propulsive  capture 
maneuver.  While  aerocapture  greatly  reduces  the  propellant  mass  required  for  orbit 
insertion,  the  large  heat  loads  generated  in  the  atmospheric  portion  of  the  trajectory 
require  the  addition  of  a  possibly  massive  thermal  protection  system  (TPS).  Both  low 
thrust  and  aerocapture  trajectories  are  highly  non-linear  with  no  closed-form  solutions  for 
their  optimal  open-loop  control  histories.  This  thesis  addresses  the  feasibility  of 
simultaneously  optimizing  a  low  thrust  interplanetary  trajectory  with  a  terminal 
aerocapture  maneuver  to  maximize  the  final  scientific  payload  mass.  In  addition,  the 
characteristics  of  optimal  low  thrust  and  aerocapture  trajectories  are  explored  with 
particular  emphasis  on  techniques  for  the  verification  of  the  optimality  of  the  obtained 
solutions. 
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D.  LOW  THRUST  MODELING 


A.  COORDINATE  SYSTEM 

Polar  coordinates  were  chosen  for  the  low  thrust  portion  of  the  trajectory  design 
optimization.  With  the  exception  of  the  planets  Mercury  and  Pluto,  the  orbits  of  the 
remaining  planets  all  lie  within  3.4°  inclination  of  the  solar  ecliptic.  For  the  purposes  of 
this  thesis,  the  small  angle  approximation  has  been  used  to  assume  that  all  planets’  orbits 
lie  in  the  solar  ecliptic.  Figure  1  depicts  the  coordinate  system  used  for  the  low  thrust 
portion. 
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Figure  1:  Low  Thrust  Coordinate  System 
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B.  EQUATIONS  OF  MOTION 


The  state  vector  of  a  vehicle  conducting  a  generic  thrusting  orbit  transfer  can  be 
given  as  x  =  [r  ,0 ,  vr,  vt,  m ] 1  where, 

r  =  radial  position  from  central  body 
0  =  transfer  angle 

v,  =  radial  velocity  component 

vt  =  transverse  velocity  component 

m  =  vehicle  mass 

The  equations  of  motion  for  a  constant  thrust  orbit  transfer  are  given  in  [Ref.  1] 
and  can  be  modified  for  non-continuous  thrust  as  follows: 


(1) 


_v, 

dt  r 


(2) 


dvr  _v,2  (i  T  sinri 
dt  r  r 1  m 


(3) 


dvt  _  vrvt  T  cost) 
dt  r  m 


(4) 


dm  _  T 
dt  ve 

where 


(5) 


(l  =  gravitational  constant  of  central  body 

T  =  thrust  magnitude 

T]  =  thrust  direction  angle 

ve  =  exhaust  velocity 
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c. 


GRAVITATIONAL  MODEL 


Only  the  gravitational  effects  of  the  central  body  are  considered  in  this  work. 
Higher  fidelity  models  include  the  gravitational  effects  of  other  celestial  bodies  thus 
enabling  the  discovery  of  gravity  assist  trajectories.  For  the  case  of  general  planet-to- 
planet  orbit  transfer,  the  vehicle  is  far  from  other  gravitational  bodies  for  the  bulk  of  the 
trajectory.  By  assuming  that  the  trajectory  begins  and  ends  at  the  boundary  of  the  sphere 
of  influence  (SOI)  with  an  escape  velocity  greater  than  or  equal  to  zero,  the  gravitational 
effects  of  the  origin  and  target  planet  can  be  neglected  completely. 


For  the  inner,  less  massive  planets,  the  SOI  are  reasonably  small  when  compared 
to  the  orbit  radii.  Eqn.  (6)  [Ref.  2]  approximates  SOI  radius  as  a  function  of  the  mean 
orbit  radius  and  the  ratio  of  masses  between  the  gravitating  bodies. 


m 


planet 


m.. 


'  sun— planet 


(6) 


Note  that  rS0I  approximately  scales  linearly  with  the  distance  to  the  sun.  Despite  the  fact 

that  the  SOI  radii  of  the  other  planets  increase  with  distance  from  the  sun,  this  increase  is 
roughly  proportional  to  the  increase  in  trajectory  length.  Considering  this,  it  is 
reasonable  to  assume  that  the  origin  and  target  planet  can  be  thought  of  as  non-gravitating 
point  masses,  reducing  the  problem  to  a  two-body  problem.  Table  1  gives  the 
approximate  SOI  radius  for  selected  planets. 
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Planet 

SOI  Radius  (km) 

SOI  Radius  (AU) 

Mercury 

1.13x10  s 

7.55xl0~4 

Venus 

6.17x10s 

4.12x10  3 

Earth 

9.24x10s 

6.18xKT3 

Mars 

5.74x10s 

3.84xl0“3 

Jupiter 

4.83xl07 

0.32 

Neptune 

8.67xl07 

0.58 

Table  1:  SOI  Radius  for  Selected  Planets  (adapted  from  [Ref.  2]) 


D.  LAUNCH  VEHICLE  MODEL 

A  launch  vehicle’s  performance  is  generally  measured  by  the  amount  of  mass  it 
can  place  in  a  given  orbit  around  the  origin  planet.  For  interplanetary  trajectories,  we  are 
concerned  with  the  trade-off  between  launch  mass  and  the  square  of  the  escape  velocity 
(know  cryptically  as  C3).  [Ref.  2].  The  launch  performance  may  be  fitted  to  a  curve  to  be 
used  to  estimate  launch  mass  from  Q  and  vice  versa.  For  this  work,  all  trajectories  will 
begin  at  the  zero  C3  point  (i.e.  maximum  launch  mass  with  escape  energy,  equivalent  to  a 
parabolic  escape). 


E.  PROPULSION  MODEL 

The  low  thrust  propulsion  system  used  for  this  research  is  modeled  after  the 
NASA  Solar  electric  propulsion  Technology  Applications  Readiness  project  (NSTAR). 
This  engine,  employed  in  the  Deep  Space  I  (DS1)  mission  [Ref.  4],  uses  a  single  xenon 
ion  propulsion  system.  This  engine  is  30  cm  in  diameter  with  a  mass  of  approximately  8 
kg.  The  engine  delivers  approximately  92  mN  of  thrust  with  a  specific  impulse  of  3300 
seconds.  To  achieve  higher  thrust  levels,  multiple  engines  can  be  used.  The  engine’s 
specific  impulse  varies  with  the  output  thrust  which  is  itself  a  function  of  the  electrical 

power  delivered  to  the  engine.  This  model  was  simplified  to  assume  a  constant  specific 
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impulse  over  the  entire  thrust  range.  Moreover,  no  limitations  were  placed  on  the 
available  power  provided  to  the  engine.  This  assumption  is  valid  provided  a  “near¬ 
limitless”  source  of  power  is  available  such  as  Nuclear-Electric  Propulsion  (NEP). 

The  propulsion  model  can  be  modified  to  account  for  “limited”  power  sources 
such  as  Solar- Electric  Propulsion  (SEP)  by  constraining  the  power  available  to  the  engine 
using 


T)  PA 

1  \  solar1  0 1  array 


1  avail 


(7) 


where 

P avail  =  electrical  power  available  to  the  engines 

Po  =  reference  solar  incident  power  (power  incident  at  1  AU)  in  W/nr 
Acn-ay  =  solar  array  area  in  m2 

T|  Soiar  =  power  converting  efficiency  of  the  solar  array 

rau  =  vehicle  distance  from  the  sun  (in  AU) 

The  power  available  can  then  be  used  to  constrain  the  maximum  thrust  available  at  a 
given  distance  from  the  sun. 


F.  NON -DIMENSIONALIZATION 

Non-dimensionalization  (or  scaling)  is  the  reformulation  of  the  problem  such  that 
all  the  optimizable  parameters  as  well  as  the  cost  function  assume  values  close  to  unity 
over  most  of  the  domain  of  interest.  When  not  scaled  properly,  the  problem  may  become 
“ill-conditioned”  and  the  effects  of  numerical  artifacts  as  common  as  round-off  error  can 
begin  to  impact  the  solution.  Very  poorly  scaled  problems  may  even  behave  as  though 
there  are  singularities  when  in  fact  there  are  none.  Problems  that  have  not  been  scaled 
are  susceptible  to  various  ailments  ranging  from  slow  convergence  to  complete  failure  to 
converge.  Betts  [Ref.  5]  gives  many  tips  for  the  scaling  of  problems  and  points  out  that 
even  when  a  problem  is  well  scaled  at  one  point,  it  may  be  scaled  poorly  at  another  point. 
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A  simple  non-dimensionalization  method  is  to  simply  divide  a  state  variable  x  by 
a  representative  reference  value  for  that  variable  [Ref.  6].  The  low  thrust  problem  may 
be  normalized  in  this  way  using  a  modified  version  of  the  approach  implemented  by 
Bryson  [Ref.  1],  Let  a  distance  unit  be  defined  as  1  AU  and  let  a  velocity  unit  be  defined 
as  the  circular  velocity  of  the  earth.  Thus  we  have 


where 

Udist  =  distance  unit  (1  AU) 

UVei  =  velocity  unit  (circular  velocity  at  1  AU) 

Using  this  convention,  a  given  radius  rx  could  be  non-dimensionalized  as  follows 


(10) 


where  an  over-bar  has  been  used  to  denote  a  non-dimensional  variable.  Continuing,  the 
natural  unit  of  time  is  simply 

U  ^y_dist_  (ii) 

^  time  jj  V A  A  / 

U  vel 

where 


Utime  =  time  unit 

Note  that  for  the  above  choice  of  distance  and  velocity  units  there  are  2t  time 
units  per  orbit.  Now  choosing  the  initial  vehicle  mass  as  our  mass  unit  we  have 

Umass=m 0  (I2) 

where 


Um  =  mass  unit 
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which  completes  the  set  of  fundamental  units.  These  particular  choices  of  units  are 
known  collectively  as  canonical  units.  [Ref.  7]  One  advantage  of  using  canonical  units 
for  astrodynamic  problems  is  that  the  gravitational  parameter  p  becomes  one, 

simplifying  many  equations.  A  second  advantage  is  that  if  only  the  above  units  (or 
combinations  of  the  above  units)  are  used  to  dimensionalize  all  the  states,  then  the 
equations  of  motion  take  on  exactly  the  same  form  as  their  dimensional  counterparts 
provided  there  are  no  constant  factors  in  any  term  which  must  be  corrected  (as  is  the  case 
with  the  gravitational  parameter  which  fortuitously  is  “corrected”  to  one).  While  it  is  true 
that  the  combinations  of  the  above  relations  can  be  used  to  non-dimensionalize  any 
parameter,  there  is  no  guarantee  that  the  result  will  also  be  well  scaled.  In  fact,  the  low 
thrust  problem  illustrates  this  fact  nicely  for  particularly  low  thrusts. 

Consider  a  600  kg  vehicle  with  one  NSTAR  engine  producing  92  mN  of  thrust. 
Using  the  units  above,  it  would  seem  reasonable  to  assume  that 


U 


thrust 


_  ^  mas  s^  dist 

U,J 


(13) 


Substituting  the  values  of  the  canonical  we  get 


U 


thrust 


( 600kg )  ( 1 .496  x  1 01 1  m ) 

- - — - - - -  =  3.9N 

(5.023x10 6)~ 


(14) 


Thus  for  this  to  be  a  useful  unit  for  non-dimensionalizing  thrust,  we  require  our  vehicle’s 
thrust  output  to  be  on  the  order  of  4  N.  For  realistic  low  thrust  systems,  this  value  of 
thrust  unit  may  be  too  high.  A  better  choice  would  be  to  define  a  thrust  unit  as  the 
maximum  thrust  produced  by  the  engine.  This  guarantees  that  the  normalized  thrust 
magnitude  will  fall  between  zero  and  one  but  requires  the  equations  of  motion  to  be 
modified  slightly  as  follows: 


(15) 


d&  _v, 

dt  7 


(16) 
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(17) 


where 


dvr  _  v,2  1  +  zT  sinri 

dt  7  7“  m 

dv,  _  _  vrv,  zT  cosrf 
dt  7  in 

dm  _  zT 
dt  v 

e 


^  thms/d  time 

dist  ^d Vel 


(18) 

(19) 


(20) 


Clearly  multiplying  all  the  non-dimensional  thrust  terms  in  the  equations  of 
motion  by  this  factor  effectively  accomplishes  the  same  objective  as  normalizing  thrust 
with  purely  canonical  units.  However,  the  difference  is  that  domain  of  the  thrust  control 
is  now  well-scaled,  and  the  conversions  are  taking  place  internal  to  the  equations  where 
proper  scaling  is  not  as  important. 
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HI.  AEROCAPTURE  MODELING 


A.  COORDINATE  SYSTEMS 

The  coordinate  systems  chosen  for  aerocapture  are  the  same  as  those  used  in 
existing  developed  tools  such  as  ACAPS  [Refs.  8,9]  such  that  data  could  more  easily  be 
shared  between  programs. 

1.  Position  Vector 

The  vehicle’s  position  vector  is  defined  in  the  planet  centered  fixed  (PCF) 
reference  frame.  This  coordinate  system  is  similar  to  a  planet  centered  inertial  (PCI) 
coordinate  system  except  that  the  primary  axis  is  aligned  and  rotates  with  a  meridian  line 
on  the  planet’s  surface  [Ref  7].  The  angular  velocity  of  the  PCF  frame  with  respect  to  the 
PCI  frame  is  denoted  by  £2.  As  shown  in  Figure  2,  the  position  vector  is  defined  as 

r  =  [r.G,^]'  where 

r  =  distance  between  origin  and  vehicle  center  of  mass 
0  =  longitude  as  measured  positive  in  an  easterly  direction 

(f>  =  latitude  measured  positive  up  from  the  equator 
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Figure  2:  Aerocapture  Position  Coordinate  System 


Because  the  targeting  of  orbits  with  specific  right  ascensions  is  beyond  the  scope 
of  this  thesis  and  the  equations  of  motion  are  invariant  in  longitude,  it  is  convenient  to 
define  the  longitude  at  atmospheric  entry  to  be  zero. 


2.  Velocity  Vector 

The  velocity  vector  is  measured  in  polar  coordinates  relative  to  the  REN  (Radius - 
East-North)  frame.  As  shown  in  Figure  3,  the  REN  frame  is  centered  on  the  vehicle  with 
the  primary  axis  aligned  with  the  extended  radius  vector.  The  tertiary  axis  is  orthogonal 
and  aligned  toward  the  northern  pole  of  the  spherical  coordinate  system  and  the 
intermediate  axis  is  mutually  orthogonal  and  oriented  easterly,  completing  the  right- 
handed  triad. 
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Radial 


Figure  3:  Aerocapture  Velocity  Coordinate  System  (REN  Frame) 

Within  this  frame,  the  velocity  vector  is  represented  by  a  magnitude  and  two 
angles,  \PCF  =  [v,\|/,y]T  where 

v  =  vehicle  speed 

\)/  =  heading  angle  as  measured  counter-clockwise  from  easterly 

y  =  flight  path  angle  as  measures  positive  up  from  E-N  plane 

It  is  important  to  note  that  the  velocity  vector  represented  in  this  frame  is  non- 
inertial  unless  Q,  =  0  in  which  case  PCF  =  PCI.  For  reasons  that  will  become  apparent 
later,  it  is  sometimes  convenient  to  express  the  vehicle’s  inertial  velocity  in  the  local 

REN  frame.  This  will  be  denoted  by  VPCT  =[V,vP,r]1  where  capital  Greek  letters  have 
replaced  their  corresponding  lower-cased  non-inertial  symbols. 
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3. 


Frenet  Frame 


The  Frenet  frame  (sometimes  referred  to  as  the  NTW  frame  [Ref  7])  is  useful  for 
resolving  vectors  into  tangent,  normal  and  bi -normal  components  with  respect  to  the 
vehicle’s  path.  At  a  given  time,  the  tangent  component  is  tangent  to  the  vehicle  trajectory 
(alternately  parallel  to  the  velocity  vector).  The  normal  component  lies  in  the  plane  of 
the  orbit  but  perpendicular  to  the  tangent  component  and  the  bi-normal  component  is 
normal  to  the  orbit  completing  the  right-hand  triad.  These  components  are  identified  by 
the  subscripts  s,  n,  and  w  respectively.  The  Frenet  frame  is  a  natural  choice  for  resolving 
aerodynamic  forces  on  a  body  as  drag  always  acts  in  the  negative  direction  tangent  to  the 
flight  path  and  the  lift  vector  will  always  lie  in  the  normal/bi-normal  plane.  [Ref  7]  That 
is  to  say  that  for  a  rotating  atmosphere  (fixed  to  the  PCF  frame)  the  Frenet  axes  are 
equivalent  to  wind  axes. 


B.  EQUATIONS  OF  MOTION 

For  the  generation  of  optimal  trajectories,  a  three  degree  of  freedom  (DOF)  model 
is  sufficient  to  describe  the  spacecraft  dynamics.  The  state  of  the  vehicle  is  represented 

by  the  7-dimensional  vector  x  =  [r,0  ,<[>  ,v  »/  ,y]'  consisting  of  radial  position,  longitude, 

latitude,  speed,  heading  angle,  and  flight  path  angle.  The  equations  of  motion  for  a  non¬ 
rotating  atmosphere  with  a  generic  gravity  model  are  given  by  [Ref.  10]  as 


dr 

—  =  vsiny 
dt 

(21) 

<r/0  v cosy  cost)/ 
dt  rcoscf) 

(22) 

d§  vcosysint)/ 

dt  r 

(23) 

dv 

~ -  =  as  +  Ss 
dt 

(24) 

aw  +  8w  _  ^_cos y  cos\j/  tan  (f> 
vcosy  r 

(25) 
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(26) 


dy  _a„+g„  [  vcosy 
dt  v  r 


where  in  the  above  equations  as ,  an ,  and  aw  are  the  external  accelerations  resolved  in 
the  Frenet  frame. 


1.  Rotational  Effects 

For  the  higher- fidelity  case  of  a  rotating  atmosphere,  Eqs.(9)  through  (11)  must  be 
modified  to  account  for  non-inertial  centrifugal  and  Coriolis  forces. 


dv 

dt 


=  as  +  gs+cfv 


d\tf  a  +  s  v  „ 

— -  =  — - —  —  cos  y  cost)/  tan(f>  +  cf  +  co 

dt  vcosy  r  ' 


dy  a„+g„  vcosy  . 

—  =  + - -  +  cf  +CO 

dt  v  r 


(27) 

(28) 

(29) 


where  the  cf  and  co  terms  represent  centrifugal  and  Corilolis  contributions  respectively  as 
given  by 


cf,  =  Q2rcost|)  (siny  cos<|)  -cosy  sintf)  sintj/ ) 


cL  =  - 


Qrr 

vcosy 


•sin  (f)  cos(]>  cost)/ 


co  =2f2(tany  cos(|)sint)/  —  sin  cf) ) 


Q?r 


cf  = - cost])  (cosy  cost])  +  siny  sintf)  sint)/  ) 


co  =  2f2cost|)  cost)/ 


(30) 

(31) 

(32) 

(33) 

(34) 


and  Q,  is  the  planet’s  angular  rotation  rate  about  its  pole. 
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2.  Aerodynamic  Forces 

The  aerodynamic  forces  of  lift  and  drag  are  given  [Ref.  1 1]  as 


L=q„SCL(a)  (35) 

D  =  qdSCD(a)  (36) 

where  Cl  and  Cd  are  the  lift  and  drag  coefficients  and  are  in  general  a  function  of  angle 
of  attack  (a  )  and  S  is  a  reference  area  associated  with  the  aerodynamic  coefficients. 
The  angle  of  attack  is  defined  as  the  angle  between  vehicle’s  body  axis  and  the  velocity 
vector  as  shown  in  Figure  4, 


Radial 


Figure  4:  Definition  of  Angle  of  Attack  and  Flight  Path  Angle 


The  dynamic  pressure,  qc\  is  given  by 
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(37) 


Qd  =^p(r)v2 

where  p  is  the  atmospheric  density  which  is  generally  a  function  of  altitude  as  given  by 
the  atmospheric  model. 

3.  External  Accelerations 

The  terms  ,  an ,  and  aw  represent  the  external  accelerations  resolved  in  the 
tangential,  normal  and  bi-  normal  directions.  These  accelerations  are  the  result  of  the 
aerodynamic  forces  and  vary  with  bank  angle,  8  which  is  defined  as  a  body  rotation  about 
the  vehicle’s  velocity  vector  as  show  in  Figure  5. 


Velocity 
(or  tangent) 


Bi-normal 

Figure  5:  Definition  of  Bank  Angle 
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These  accelerations  can  be  expressed  as  function 

of  the  controls,  mass  and 

aerodynamic  forces  as 

D 

(38) 

as  = - 

m 

Lc  os  8 

(39) 

a»  = 

m 

L  sinS 

(40) 

a  = - 

W 

m 

Note  that  the  total  acceleration  on  the  vehicle  can 

be  determined  from  these 

equations  by  using 

\\a\\  =  \las2  +  an2  +  aw2 

(41) 

and  the  total  g-load  becomes 

g  -  load  =  U 

(42) 

So 

C.  GRAVITATIONAL  MODEL 

For  a  simple  inverse- square  gravity  model,  the  gravitational  accelerations  can  be 

converted  from  their  more  familiar  spherical  form 

ft 

g=gr  = - V 

r 

(43) 

go  =0 

(44) 

g*  =0 

(45) 

to  the  Frenet  frame  using  Eqn.  (215)  becoming 

8n  =-gcosy 

(46) 

gs  =  ~g  sinY 

(47) 

gw  =o 
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(48) 

Increased  fidelity  gravitational  models,  such  as  those  including  zonal  harmonics 
can  be  incorporated  in  a  similar  manner  [Ref.  9]  but  is  beyond  the  scope  of  this  work. 

D.  ATMOSPHERIC  MODEL 

Anderson  shows  [Ref.  11]  that  if  an  atmosphere  is  assumed  to  be  isothermal  with 
constant  temperature  equal  to  the  mean  temperature,  then  atmospheric  density  varies 
exponentially  with  altitude.  That  is, 

(ft-jo) 

P(r)  =  Poe  Hp  (49) 

where 

p o  =  reference  density 

Hp  =  scale  height 

h  =  altitude  (  r-rpianet) 

ho  =  reference  altitude 

The  atmosphere  is  assumed  to  be  fixed  to  the  planet  and  rotates  with  the  PCF 
coordinate  system  about  the  KPCF  axis.  The  following  table  provides  the  above  data  for 
Mars  and  Neptune. 


Parameter 

Mars  Value 

Neptune  Value 

Po 

4.7x10  ^kg/m3 

2.348x10 r3kg/m3 

HP 

l.OxlO4  m 

5.331xl04  m 

h0 

4.9xl04  m 

0  m 

Atmosphere  Limit 

12  5  km 

8  00  km 

Table  1:  Atmospheric  Parameters  for  Mars  and  Neptune 
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E. 


HEATING  MODEL 


The  thermal  protection  system  (TPS)  for  a  re-entry  vehicle  is  typically  sized  to 
sustain  both  the  heating  rate  and  heat  load  (integral  of  heating  rate)  encountered  by  the 
spacecraft  [Refs.  12,  13]  These  quantities  take  a  maximum  value  at  the  stagnation  point, 
the  point  on  the  surface  of  the  vehicle  where  the  freestream  velocity  is  brought  locally  to 
zero.  The  heating  rate  at  this  point  can  be  approximated  by  the  Chapman  equation 


q  = 


(50) 


where 


rn  =  vehicle  nose  radius 


C  =  stagnation  point  heating  coefficient 

and  N  and  M  are  constants  for  a  given  atmosphere  associated  with  density  and  speed 
respectively.  [Ref.  14] 


Parameter 

Mars  [Ref.  15  ] 

Neptune  [Ref.  17] 

C 

3.55xl05 

7.9x10 5 

N 

0.5 

0.5 

M 

3.15 

3.0 

Table  2:  Heating  Rate  Parameters  for  Mars  and  Neptune 


F.  VEHICLE  MODEL 

The  vehicle  data  used  was  obtained  from  [Ref.  9]  for  consistencies  with  that 

author’s  original  work  in  aerocapture  and  is  based  loosely  on  the  Mars  Pathfinder  shape. 

Due  to  the  high  cost  to  flight  qualify  new  space  flight  components,  most  re-entry  heat- 

shields  have  deviated  very  little  since  the  early  Viking  designs.  In  light  of  this  fact,  the 

aerodynamic  data  below  will  likely  be  representative  of  possible  future  aerocapture 

vehicles.  The  aerodynamic  coefficients  listed  below  apply  to  Mars  trajectories;  lift  and 
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drag  coefficients  at  other  planets  should  vary  primarily  due  to  variations  in  each  planets’ 
atmospheric  composition.  Unfortunately  aerodynamic  data  for  a  Neptune  aerocapture 
mission  is  not  available;  therefore  the  values  for  Mars  will  be  used  with  the 
understanding  that  these  coefficients  do  not  agree  with  the  vehicle  design. 

Following  the  atmospheric  pass,  the  vehicle  will  require  a  small  impulsive  AV  at 
apoapsis  to  place  the  vehicle  in  its  final  circular  orbit.  To  accomplish  this,  a  small 
conventional  chemical  engine  using  a  bi-propellant  comprised  of  Hydrazine  (N2H4)  fuel 
and  Nitrogen  Tetroxide  (N2O4)  as  oxidizer  is  employed.  This  combination  yields  a 
specific  impulse  of  approximately  330  seconds  [Ref.  13]  and  provides  sufficient  thrust 
such  that  the  AV  may  be  considered  impulsive. 


Parameter 

Symbol 

Value 

Vehicle  Mass 

m 

568.5  kg 

Nose  Radius* 

rn 

1  m 

Coefficient  of  Lift 

CL 

0.3024 

Coefficient  of  Drag 

Cd 

1.6800 

Reference  Area 

S 

5.52  nr 

Engine  Specific  Impulse 

Isp 

330  sec 

Engine  Exhaust  Velocity 

Ve 

3234  m/s 

Table  3:  Aerocapture  Vehicle  Parameters 

G.  NON -DIMENSIONALIZATION 


The  aerocapture  equations  of  motion  have  been  non-dimensionalized  in  the  same 
manner  used  by  [Ref.  9]  in  his  original  work  in  aerocapture  and  are  repeated  here  for 
completeness. 


*  Note:  This  value  is  not  the  actual  measure  of  the  nose  radius;  rather  this  value  is  chosen  to  be 
consistent  with  the  stagnation  point  heating  coefficient  C  in  Table  3 
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_  r 

r 


U 


dist 


where  U dist  =  r  ,t ,  the  planet’s  surface  radius. 


U 


vel 


where  Uvel  =  __  — — ,  the  circular  velocity  at  the  planet’s  surface. 


planet 


(51) 


(52) 


where  Utime  =^-^l  and 


U, 


vel 


^ time 


(53) 


_  m 
m  —  ■ 


U„ 


(54) 


where  Uinass  =m0 ,  the  vehicle’s  initial  mass.  Noting  that  the  angles  measured  in  radians 

are  already  conveniently  non-dimensional  and  well  scaled,  Eqns.  (21)-(23)  may  be 
presented  in  their  non-dimensional  form 


dr 

dt 


■  v  siny 


dQ  _v  cosy  cost)/ 
dt  7  c  os(f> 

d§  _  v  cosy  sint \f 
dt  7 


(55) 

(56) 

(57) 


The  aerodynamic  forces  may  be  non-dimensionalize  by  first  assuming  that 

^  P 


U 


(58) 


density 
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where  U demity  =  p0 ,  the  scale  density  and  using  a  non-dimensional  reference  area  defined 
as 


5  = 


Noting  that  the  dynamic  pressure  non-dimensionalizes  as 


f  U  A 

mass 

i  ^ density^ dist  j 


the  aerodynamic  forces  may  be  non-dimensionalized  as 

L=qdSCL{a) 

D  =  qdSCD  (oc ) 

The  external  accelerations  may  now  be  written  as 

D 

as  =-  — 

m 

_  LcosS 

an  = - _ 

m 

_  LsinS 

a  = - 

W  — 

m 

The  gravitational  acceleration  becomes 


where  Ugrav  =—^LT. 

^ planet 

Eqns.  (27)-(29)  become 


8  =  gr 


U 


grav 


dv  _  _ 

~=  =  as+gs+cfv 
dt 


(59) 

(60) 

(61) 

(62) 

(63) 

(64) 

(65) 

(66) 


(67) 
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(68) 


chit  au.  +  gw  v  ,  _ 

—hr =~z - cosy  cost)/  tm^+cf  +co 

dt  vcosy  r 


dy 

dt 


a,+g„  vcosy  _ 

+  — +  cf  +  CO 
v  r 


(69) 


where  the  centrifugal  and  Coriolis  forces  have  been  non-dimensionalized  using 

n=nu[ime  (70) 

such  that 


Cl2r  cos(|)  (siny  cos(|)  -  cosy  sintf)  sint)/ ) 

(71) 

H2r 

cf  = - sin(f>  costf)  cost)/ 

(72) 

v  cos  y 

cbv  =  2£2(tany  cos(|)sin\|/  —  sin  (f> ) 

(73) 

D.2t 

cos (f>  (cosy  cos ([)+  siny  sin(|)  sint)/ ) 

(74) 

v 
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IV.  SOLVING  OPTIMAL  CONTROL  PROBLEMS 


A.  PRELIMINARIES 

An  optimal  control  problem  is  the  task  of  solving  for  the  state  and  control 
histories  of  a  system  subject  to  constraints  while  minimizing  (or  maximizing)  some 
performance  index.  Put  mathematically  (as  presented  by  [Ref.  18]  and  repeated  here  for 
completeness),  given  a  system  with  dynamic  constraints  such  that 

x(x)=f(x,u,x;p)  (75) 

where 

x  =  vector  of  states  that  completely  describe  the  systems  at  any  x 
u  =  vector  of  controls 

x  =  independent  variable  (usually  but  not  necessarily  time) 
p  =  vector  of  static  parameters 
f  =  vector  of  dynamic  equations 
subject  to  additional  path  constraints 

g/<g(x,u,x)<g„  (76) 

g  =  vector  path  constraint  equations 

gi  =  vector  of  lower  path  bounds 
gu  =  vector  of  upper  path  bounds 
and  boundary  conditions  (or  point  constraints) 

e,  <e(x(x0),x(x/),x0,x/)<e„  ill) 

e  =  vector  path  event  condition  equations 
ei  =  vector  of  lower  event  condition  bounds 

eu  =  vector  of  upper  event  condition  bounds 
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as  well  as  bounds  on  states  and  controls 

x;<x(x)<xu  (78) 

u,<u(x)<u„  (79) 

in  order  to  minimize  a  performance  index  of  form 

x/ 

J  (x(»),u(»),T0,T/)  =  £'(x(x0),x(x  f  ),x0,x  y)  +  Jf(x(x),u(x),x  )d%  (80) 

Xo 

E  =  scalar  cost  function  evaluated  at  the  boundary  times  (event  cost) 

F  =  scalar  cost  function  evaluated  over  the  entire  time  history  (integral  cost) 

When  Fe  0,  the  cost  function  is  said  to  be  in  Mayer  form.  When  Ee0 ,  the  cost 
function  is  said  to  be  in  Lagrange  form.  When  E,F  <t0  the  cost  function  is  said  to  be  in 

Bolza  form. 

A  mathematical  construct  known  as  the  Hamiltonian  may  be  created  by  adjoining 
a  vector  of  Lagrange  multipliers  to  the  dynamic  constraints  and  adding  the  integral  cost 

H  (x,u,A,,x;p)  =  F(x(x),u(x),x;p)+XT*f  (x(x),u(x),x;p)  (81) 

(Note:  the  Lagrange  multipliers  associated  with  the  dynamics  are  also  referred  to  as 
“costates”).  The  Pontryagin  Minimum  Principle  (PMP)  [Ref.  1]  states  that  the  optimal 
control  history  satisfies 

u  =  argmin  H  subject  toucU  (82) 

where  u  is  the  optimal  control  and  U  is  the  domain  of  u. 

This  should  be  recognized  as  a  static  optimization  problem  at  each  instant  of  time 
for  which  the  Hamiltonian  itself  is  now  the  performance  index  minimized.  We  can  now 
form  an  additional  construct  known  as  the  augmented  Hamiltonian  (or  Lagrangian  of  the 
Hamiltonian)  defined  as 

H 1  (p,?i,x,u,x;p)  =H (x,u,?i,x;p)+pT-g  (x(x),u(x),x;p)  (83) 
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The  Karush-Kuhn-Tucker  conditions  state  that  when  evaluated  at  the  extremal 


controls,  the  following  are  true 


du 

MTg  =  0 


<0 

gi  =  8 

II  IV 
o  o 

8  =  8U 

&0 

VI 

so 

VI 

^0 

any 

3 

^0 

II 

&0 

(84) 

(85) 


(86) 


where  Eqns.  (84)  and  (85)  are  referred  to  as  the  gradient  normality  and  complementary 
conditions  respectively.  Eqn.(86)  will  be  useful  later  for  verifying  the  switching  structure 
of  the  controls. 


B.  SOLUTION  METHODS 

Methods  for  solving  optimal  control  problems  can  generally  be  separated  into  two 
groups,  indirect  and  direct  methods  [Ref.  5  p.  85]. 

Indirect  methods  tend  to  produce  greater  accuracy  and  faster  solution  times. 
However,  the  problem  formulation  for  an  indirect  method  is  considerably  more  complex 
as  the  user  must  provide  additional  information  such  as  the  equations  for  the  costate 
dynamics  as  well  as  the  gradient  of  the  Hamiltonian  with  respect  to  the  controls. 
Moreover,  one  must  provide  a  reasonably  accurate  guess  for  the  controls,  states  and 
costates.  Unfortunately,  the  abstract  nature  of  the  costates  makes  this  guess-work 
something  of  an  art. 

By  contrast,  direct  methods  tend  to  be  much  more  robust  and  can  converge  to  an 
optimal  solution  even  when  seeded  with  a  poor  cr  infeasible  guess.  Direct  methods 
require  less  prior  work  on  the  part  of  the  user  as  the  user  need  only  supply  the  dynamic, 
path  and  event  constraint  functions  as  well  as  the  cost  function.  Direct  methods  work  on 
the  premise  that  a  continuous  problem  may  be  approximated  through  careful 
discretization  as  a  large  number  of  point  constraints,  thus  reducing  the  problem  to  a 
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single  large  Nonlinear  Programming  problem  (NLP).  Unfortunately,  to  increase  the 
accuracy  of  the  solution,  greater  resolution  is  required  leading  to  larger  problems  that  can 
be  very  slow  to  converge. 


C.  DISCRETIZATION 

Approximation  theory  tells  us  that  the  worst  (with  regard  to  minimizing  the  L 
error  norm)  discretization  scheme  is  to  sample  a  function  at  equal  intervals  [Ref.  6]. 
Other  methods  of  discretization  include  Hermite-  Simpson  and  Sine  methods  [Ref.  20].  In 
fact,  discretizing  a  function  at  the  Legendre-Gauss-Labatto  (LGL)  points  minimizes  the 
Lo  error  norm  for  a  given  number  of  nodes.  LGL  points  are  characteristic  in  that  they 
have  a  higher  density  distribution  at  the  end  points  of  a  function  and  becoming  sparser  in 
the  interior  regions.  By  discretizing  the  problem  in  this  manner,  a  direct  solution  may  be 
obtained  with  either  more  accuracy  or  less  computation  time. 


D.  DIDO 

DIDO  is  an  application  package  [Ref.  18]  for  solving  dynamic  optimization 
problems  in  a  friendly,  easy-to-use  MATLAB  environment.  DIDO  employs  a  powerful 
direct  Legendre  pseudospectral  method  that  exploits  the  sparsity  pattern  of  the  discrete 
Jacobian  by  way  of  the  NLP  solver  SNOPT”  [Ref.  21]  and  runs  in  both  UNIX  and  PC 
environments. 

Problem  formulation  in  DIDO  is  quite  simple,  with  the  user  creating  a  set 
MATLAB  functions  to  evaluate  the  problem  dynamics,  cost,  path  constraints  and  event 
conditions.  These  functions  are  tied  together  by  an  additional  script  which  defines  the 
upper  and  lower  bounds  for  the  states,  controls,  path  constraints  and  event  conditions  as 
well  as  the  initial  guess.  This  guess  need  not  be  feasible;  in  most  cases  simply  providing 
the  estimated  initial  and  final  state  value  is  sufficient.  The  simplicity  of  DIDO  contributes 
to  the  rapid  prototyping  of  solutions.  An  experienced  DIDO  user  can  generate  optimal 
trajectories  in  days  that  previously  would  have  required  weeks  or  months. 
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The  user  defines  the  number  of  nodes  (discretization  points)  for  the  problem  and 
is  afforded  the  option  of  “cold  starting”  or  “bootstrapping”  from  a  previous  solution. 


E.  COVECTOR  MAPPING  THEOREM 

Another  advantage  of  the  Legendre  pseudospectral  method  is  that  is  provides 
accurate  costate  and  covector  histories  despite  the  fact  that  the  adjoint  equation  are  not 
supplied.  [Refs.  22,  23]  This  power  comes  from  the  realization  of  the  Covector  Mapping 
Theorem  (CMT).  [Ref.  21]  provides  the  following  explanation  of  the  CMT: 


The  CMT  may  be  articulated  as  follows:  Given  an  optimal  control 
problem,  P ,  let  PN  denote  its  Legendre  pseudospectral  approximation 
where  N  is  the  order  of  the  Legendre  polynomial  used  in  the 
approximation.  Let  P'  denote  the  boundary  value  problem  (BVP)  arising 
from  an  application  of  the  PMP  to  problem  P ,  and  PNX  denote  the 
generalized  root -finding  problem  obtained  by  applying  the  Karush-Kuhn- 
Tucker  (KKT)  conditions  to  problem  PN  .  The  CMT  asserts  that  if  P'  is 
discretized  by  a  Legendre  pseudospectral  method  (i.e.  an  indirect  method), 
then  there  exists  an  order-preserving  map  between  these  discretized 
covectors  and  the  KKT  multipliers  associated  with  problem  PNX .  Hence, 
from  the  KKT  multipliers,  one  can  find  covectors  by  the  direct  Legendre 
pseudospectral  method  as  if  one  solved  the  problem  by  the  indirect 
method.  This  is  the  CMT. 


The  power  of  the  CMT  is  that  it  provides  the  user  with  excellent  tools  for 
ascertaining  the  optimality  of  a  given  solution  as  we  will  see  in  the  next  section. 
Additionally,  because  the  CMT  provides  accurate  costate  information,  a  DIDO  solution 
can  be  uses  as  an  extremely  high  quality  guess  for  a  more  accurate  indirect  method. 


F.  VERIFICATION  OF  OPTIMALITY 

When  solving  optimal  control  problems  one  is  often  challenged  as  to  how  one  can 
prove  optimality.  [Ref.  21]  discusses  techniques  for  verification  of  optimal  trajectories. 
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1.  Feasibility 

Of  primal  concern  is  the  feasibility  of  a  solution.  That  is,  does  the  solution  satisfy 
the  dynamic  constrains  of  the  problem?  To  ascertain  this,  the  initial  conditions  are 
propagated  using  a  Runge-Kutta  algorithm  (the  MATLAB  command  ocle45  which 
implements  the  Dormand- Prince  pair)  using  controls  interpolated  from  the  DIDO 
solution.  If  the  DIDO  state  history  and  the  propagated  state  history  match  (within  some 
tolerance)  feasibility  is  declared. 

2.  Accuracy 

The  next  issue  addresses  the  accuracy  of  the  solution.  Assuming  the  trajectory  is 
feasible,  the  event  conditions  are  evaluated  using  either  the  DIDO  or  propagated  solution 
and  examined  to  verify  that  all  event  constrains  are  satisfied  (again,  within  a  tolerance). 

3.  Well-behaved 

Similarly,  the  solution  is  said  to  be  well  behaved  if  the  path  constraints  (including 
state  and  control  bounds)  are  obeyed  over  the  entire  trajectory.  A  solution  that  violates 
the  path  constraints  may  indeed  be  feasible;  however  it  does  not  solve  the  problem  at 
hand. 

4.  Optimality 

Finally  and  most  importantly,  does  the  solution  satisfy  the  necessary  conditions 
for  optimality?  Recalling  Eqn.  (84),  analytical  expressions  can  be  obtained  by  taking  the 
partial  of  the  augmented  Hamiltonian  with  respect  to  each  control.  In  some  cases,  these 
expressions  can  be  rearranged  into  an  explicit  control  law.  Regardless,  these  expressions 
must  hold  true  for  there  to  be  first-order  optimality.  These  expressions  can  be  verified 
through  substitution  of  the  states  and  CMT  obtained  costates  and  covectors.  In  the  case 
where  the  necessary  conditions  can  be  re-arranged  to  solve  for  a  control,  these  controls 
can  be  solved  (referred  to  as  CMT  controls)  and  compared  to  DIDO  solution  controls. 
When  the  DIDO  and  CMT  controls  are  in  agreement,  first-order  optimality  is  declared. 

A  second  method  for  determining  optimality  is  by  examination  of  the 
Hamiltonian  which  is  also  constmcted  by  DIDO  via  the  CMT.  Recalling  again  Eqn.(84), 
the  first  integral  of  the  necessary  condition  is  required  to  be  constant  for  Hamiltonians 
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that  are  not  explicit  functions  of  time.  Thus  the  Hamiltonian  can  be  plotted  and  its 
“flatness”  used  as  a  measure  of  optimality. 
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V.  OPTIMAL  LOW  THRUST  PROBLEM  FORMULATION 


A.  EVENT  CONDITIONS 

For  the  orbit  transfer  problem,  the  vehicle  must  begin  at  Earth’s  orbit  distance 

r(to)  =  rEar,h  (87) 

and  complete  the  trajectory  at  the  target  planet’s  orbital  distance 

rM  =  W t  (88) 

Assuming  circular  orbits  with  no  ephemeris,  there  are  an  infinite  number  of 
equivalent  solutions  because  the  Hamiltonian  is  invariant  with  the  vehicle’s  angular 
displacement.  Thus  it  is  convenient  to  tie  down  the  initial  angular  displacement  to  some 
value. 

e(*o)=0  (89) 


Assuming  that  the  vehicle  escapes  Earth’s  SOI  with  zero  Q,  the  vehicle’s  initial 
velocity  must  match  that  of  Earth.  Again,  assuming  circular  planet  orbits,  the  equation 
for  a  body’s  radial  and  tangential  velocities  are 

v,(»„)  =  0  (90) 


where  =  1.327x  1020 m/i  is  the  Sun’s  gravitational  parameter.  Moreover,  initial 

vehicle  mass  must  lie  between  the  maximum  capacity  for  the  launch  vehicle  and  some 
lower  design  bound. 

For  the  rendezvous  case  the  vehicle’s  final  velocity  vector  must  match  that  of  the 
target  planet  thus 

E  (tf )  =  0  (92) 
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(93) 


Sun 


'  target  planet 


For  the  case  of  a  fly-by,  the  final  velocity  components  would  be  left  unconstrained. 

B.  COST  FUNCTIONS 
1.  Minimum  Time 

To  minimize  time,  the  performance  index  can  be  constructed  in  either  Mayer  or 
Lagrange  form.  In  Mayer  form,  the  cost  function  is  simply 

J=tf  (94) 

Formulated  as  a  Lagrange  cost,  the  performance  index  becomes 

J  =  {  dr  (95) 


2.  Minimum  Fuel 

For  a  single-staged  vehicle,  minimizing  fuel  is  equivalent  to  maximizing  the  final 
vehicle  mass.  Maximizing  a  quantity  is  the  same  as  minimizing  the  negative  of  that 
quantity,  thus  the  Mayer  performance  index  for  minimum  fuel  is  simply 

J  =  —mf  (96) 


where  m/  is  the  final  vehicle  mass.  Alternately  this  performance  index  may  be  formulated 
as  a  Lagrange  cost  by  recalling  that  the  fuel  flow  rate  for  a  constant  specific  impulse 
propulsion  system  is  given  by 


m 


fuel 


(97) 


so  the  performance  index  can  also  be  expressed  as 
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(98) 


dr 


C.  CONTROLS 

The  natural  choices  of  controls  for  the  low  thrust  problem  are  thrust  magnitude 
and  thrust  angle.  However,  for  some  problems  these  choices  lead  to  additional 
difficulties.  Bounding  the  thrust  angle  as  0<T)  <  2n  can  cause  a  jump  discontinuity 
when  the  optimal  control  history  passes  through  the  boundary.  Increasing  the  bounds  to 
-2 7t  <r)  <  2k  yields  multiple  values  of  T|  that  are  equivalent,  again  possibly  contributing 
to  numerical  instability  of  the  solver  algorithm. 

To  remedy  this  problem  the  equations  of  motion  were  altered  such  that  the 
controls  are  the  radial  and  transverse  thrust  components  as  well  as  thrust  magnitude. 
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Thrust  Mag 


Figure  6:  Thrust  Cone 


From  Figure  6  clearly  the  thrust  vector  can  be  described  by  its  radial  and 
transverse  components 

Tr  -  Tmag  sinr)  (99) 

T,  =Tmag  COST) 

and  Eqns.  (3), (4)  and  (5)  can  be  rewritten  with  our  new  controls  as 

dVr  _  V  M-  ,  Tr 

dt  r  r  m 

dv<  _  vrv,  ,  T, 
dt  r  m 


(100) 

(101) 

(102) 
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(103) 


Our  control  vector  u  is  now 


dm 

dt 


where  our  controls  are  constrained  within  the  following  bounds 


-T  <T  <T 

max  r  n 

-T  <T <T 

max  t  n 

o<r  <t 

mag  ma 


(104) 


(105) 

(106) 
(107) 


D.  STATE  BOUNDS 

Each  state  may  be  limited  to  a  certain  range  of  values  over  the  trajectory,  while 
others  appear  unbounded.  As  an  example,  the  vehicle  mass  over  the  duration  of  the 
trajectory  may  not  exceed  its  initial  mass.  Assuming  that  the  vehicle’s  only  method  of 
shedding  mass  is  via  thrust,  then  the  dry  mass  of  the  spacecraft  serves  as  a  lower  bound. 
Thus  over  the  entire  problem,  mass  is  restricted  to 

mdry  ~  m  —  minitial  (108) 

which  might  take  a  non-dimensional  for  of  something  like 

.1  <m<l  (109) 

where  the  non-dimensional  value  for  the  lower  bound  depends  on  the  vehicle’s  dry  and 
initial  masses. 

Other  states  should  be  bound  for  more  practical  reasons.  For  example,  a  logical 
choice  of  bounds  for  heliocentric  radius  might  be 

0  <r<°°  (110) 
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However,  the  dynamics  [Eqns.  (2)-(4)]  become  singular  at  r  =  0,  thus  this  is  a 
bad  choice  for  lower  bound.  Similarly,  although  there  is  no  physical  reason  to  have  an 
upper  bound,  placing  a  realistic  upper  bound  of,  say,  100  reduces  the  search  space  for  the 
NLP  solver.  Thus  better  state  bounds  might  be 

.1  <  r  <  100  (111) 

where  the  bounds  have  been  expressed  in  canonical  units.  For  an  orbit  transfer  to  Mars 
(at  only  1.52  AU)  this  upper  bound  could  be  reduced  reasonably  further. 

Following  similar  reasoning,  the  remaining  states  take  non-dimensional  bounds  of 


-10  <v,,  <10 
-10  <  v,  <  10 
0 <6  <  107t 

Thus  the  state  bounds  can  be  expressed  in  vector  format  as 


"  .1  ~ 

r 

'100' 

0 

0 

10k 

-10 

< 

A 

< 

10 

-10 

v, 

10 

.1 

m 

1 

(112) 

(113) 

(114) 


(115) 


E.  PATH  CONSTRAINTS 

As  currently  bound,  the  control  space  incorrectly  takes  the  shape  of  a  rectangular 
solid.  An  additional  constraint  is  needed  to  relate  the  three  controls  and  correct  the 
control  domain. 

T* +Tt2 -T^2  =  0  (116) 

In  the  control  space  u  cz  9t3  this  path  constraint  represents  the  surface  of  a  cone. 
Because  the  interior  of  the  cone  is  not  in  the  domain  of  u,  the  problem  is  not  convex,  thus 
leading  to  additional  difficulties  with  attaining  a  rapidly  converging  solution. 
Interestingly,  this  problem  can  be  convexified  by  restating  the  path  constraint  as 
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~T  2  <t2  +j2  _j  2<q  (H7) 

x  mag  —  ±  r  '  ±  t  ±  mag  —  w  V x  x  7  / 

which  represents  a  solid  cone.  Thus  physics  must  be  violated  to  convexify  the  problem. 
Fortunately,  Pontryagin’s  Minimum  principle  [Ref.  19]  ensures  that  the  optimal  control 
lies  on  the  boundary  of  the  constraint  surface.  In  essence,  physics  is  intentionally 
violated  to  improve  the  convergence  of  the  solution  with  the  guarantee  that  the  Minimum 
Principle  will  restore  continuity  of  physical  reality. 


F.  NECESSARY  CONDITIONS 

Recall  from  chapter  IV  that  the  Hamiltonian  is  defined 

H  (x,u,A,,x;p)  =  F(x(x),u(x),x;p)+>LT»f  (x(x),u(x),x;p)  (118) 

Assuming  any  Mayer  cost  (that  is  F  e  0 ),  the  Hamiltonian  can  be  constructed  as 


H=\(vr)+\ 
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(119) 


The  augmented  Hamiltonian  H  '  is  found  by  adjoining  the  path  constraints  to  the 
above  Hamiltonian  as  follows 
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[l7, T  +|i  f  +  [i  T  +  (J.  ) 

r  lr  r  r  If  t  r  lmag  mag  r  rel  \  r  t  mag  / 


where  ,  and  u7  are  the  Lagrange  multiplies  associated  with  the  control  bounds 

1r  1  mag 

and  pre/is  the  Lagrange  multiplier  associated  with  path  constraint  relating  the  three 
controls  [Eqn.  (117)]. 

Applying  the  necessary  condition  for  optimality  by  taking  the  partial  derivative  of 
the  augmented  Hamiltonian  with  respect  to  each  control  we  obtain  the  following  three 
necessary  conditions: 

})H  X 

^T  =  — +Mt  +2^=0  (121) 

ol  m 
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dH'  A, 

—-  =  ^+^+2^=0 


(122) 


t  X 

TT  =  ~  ~  +  ~  ^  ^rePmag  =  0 

1  mag  v  e 


(123) 


These  can  be  re-arranged  to  explicitly  solve  for  the  controls 


i  (K 

Tr  =  ~zr-  — 

2  Vrel  m 


(124) 


1 

T.=-^—  —+^ 

2  M"  rel  m 


(125) 


T  =  1  [-^il+n 

"  2nJ  v,+^ 


(126) 


The  normalized  equations  of  motion  had  an  additional  factor  relating  the  thrust 
units  to  the  other  canonical  units.  Thus  in  non-dimensional  form  Eqns.  ( 1 24)-(  126) 
become 


i  f  z\ 

T,  =-tt—  -^+m-r 

2  Vrel  m 


(127) 


i  f  zk. 


2Vre, 7  m 


(128) 


T  =  1  ( -  Z^m  +  u 

2|d  v,  +  Mt- 


(129) 


where 


^  thrush  time 
^2 dist  ^  v  el 


(130) 


These  expressions  will  prove  useful  for  verification  of  the  DIDO  solutions  in  the 
next  chapter. 
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VL  OPTIMAL  LOW  THRUST  RESULTS 


A.  PROCESS 

The  code  was  first  validated  by  solving  the  well-known  minimum  time  low  thrust 
transfer  problem  as  presented  in  [Ref.  1],  With  this  success  in  hand,  the  parameters  were 
modified  to  better  represent  actual  technology.  It  was  at  this  time  that  the  problems  with 
the  original  non-dimensionaliztion  scheme  for  thrust  became  apparent.  Although 
minimum  time  trajectories  were  acceptable,  feasible  minimum  fuel  trajectories  were 
difficult  to  obtain.  The  problem  was  reformulated  using  the  improved  non- 
dimensionalization  scheme  with  better  results. 

In  general,  low  thrust  problems  took  considerably  less  computation  time  than  their 
aerocapture  counterparts  due  to  the  fewer  number  of  state  variables  and  “slower” 
dynamics.  For  this  reason,  solutions  were  generally  not  bootstrapped,  but  instead  solved 
completely  each  time  from  an  initial  guess. 


B.  MINIMUM  TIME  RENDEZVOUS 

The  minimum  time,  Mars  rendezvous  problem  was  solved  for  a  vehicle  powered 

by  six  NSTAR  ion  engines  providing  a  total  thrust  capacity  of  0.55  N.  The  initial  vehicle 

mass  was  fixed  at  the  maximum  lift  mass  (C3  of  zero)  for  the  Delta  II  7326-9.5  which 

corresponds  to  659.3  kg.  If  the  initial  mass  is  not  fixed  at  some  value,  initial  runs  showed 

that  the  solution  would  always  choose  the  smallest  initial  mass  available  with  the 

minimum  propellant  mass  necessary  to  complete  the  trajectory,  essentially  minimizing 

the  inertia  the  thruster  must  work  against.  As  discussed  in  Chapter  V,  a  zero  C3  departure 

trajectory  can  be  modeled  as  beginning  at  the  origin  planet’s  location,  matching  the 

planet’s  tangential  velocity  and  with  zero  radial  velocity.  Sixty  nodes  were  used  to  solve 

the  problem,  requiring  only  4.91  minutes  on  a  Pentium  IV  PC  (1.8  Ghz,  512  Mb  RAM). 

This  large  number  of  nodes  was  more  than  sufficient  to  capture  the  details  of  the 

trajectory  but  the  increased  number  of  nodes  was  used  to  generate  more  aesthetically 

pleasing  plots.  The  minimum  time  for  transfer  is  182.57  days  during  which  time  the 

vehicle  consumes  270.2  kg  of  propellant  for  an  arrival  mass  of  389.2  kg.  Figure  7  shows 
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the  non-dimensional  state  histories  versus  time.  The  symbols  represent  the  DIDO 
solutions  while  the  line  represents  the  propagated  solution  (using  the  DIDO  controls). 


Low  Thrust  States:  Earth  to  Mars 


Figure  7:  State  History  (Min  Time  to  Mars  Rendezvous) 


The  radius  vector  smoothly  and  asymptotically  transitions  from  1  AU  to  1.52  AU 
while  the  transfer  angle  increases  almost  linearly  to  136.6  deg.  Radial  velocity  increases 
to  a  maximum  of  0.35  canonical  units  before  symmetrically  returning  to  zero  to  effect  the 
rendezvous.  Transverse  velocity  first  slightly  increases  before  beginning  to  decay 
matching  Mars’  circular  velocity  of  0.81  canonical  units.  Mass  of  course  decreases 
linearly  due  to  the  constant  specific  impulse  and  constant  thrust  profile  as  we  shall  see 
shortly. 

The  heliocentric  transfer  orbit  is  shown  in  Figure  8  with  the  viewpoint  from  the 
solar  ecliptic  plane  north  pole.  Again,  the  circles  are  the  DIDO  trajectory  and  the  solid 
line  is  the  propagated  trajectory.  The  arrows  are  oriented  with  the  thrust  direction  and 

scaled  to  the  magnitude  of  the  thrust  vector.  There  is  a  noticeable  switch  approximately 
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half  the  distance  between  Earth  and  Mars  where  the  thrust  vector  swings  through  from 
being  predominantly  outward  to  being  predominantly  inward.  It  is  this  change  in  the 
radial  component  of  the  velocity  vector  that  begins  to  arrest  the  vehicles  radial  velocity  to 
prepare  it  for  a  zero  relative  velocity  arrival. 


Low  Thrust:  Earth  to  Mars 


Figure  8:  Heliocentric  Trajectory  (Min  Time  to  Mars  Rendezvous) 

Figure  9  shows  the  control  history  during  the  trajectory.  The  DIDO  controls  are 
plotted  as  circles  and  the  CMT  derived  controls  are  plotted  as  dots.  The  change  in  radial 
component  of  the  velocity  vector  is  readily  apparent  in  the  thrust  angle  control  history 
just  prior  to  the  100  day  mark.  As  expected,  the  DIDO  thrust  magnitude  for  the 
minimum  time  trajectory  employs  maximum  thrust  for  the  entire  trajectory.  The  DIDO 
controls  and  the  CMT  controls  are  in  excellent  agreement  for  the  thrust  anglealthough 
they  do  not  coincide  for  the  thrust  magnitude.  Despite  the  CMT  controls  telling  us 
otherwise,  the  constant  thrust  solution  is  well  known  for  the  minimum  time  problem;  thus 
it  seems  there  may  be  an  error  in  the  way  DIDO  calculates  costates  and  covectors  or  an 
artifact  specific  to  this  problem  formulation.  Figure  10  gives  the  Hamiltonian  for  the 
problem.  Note  that  the  variations  in  the  Hamiltonian  are  on  the  order  of  10'“,  that  is,  the 
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Hamiltonian  is  very  nearly  constant.  This  fact,  combined  with  the  excellent  agreement 
between  the  thrust  angle  DIDO  and  CMT  controls  suggest  that  the  first  order  necessary 
conditions  for  optimality  have  been  met. 


Control  Histories 
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Figure  9:  Control  History  (Min  Time  to  Mars  Rendezvous) 
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Hamiltonian  vs  Time 


-1.555 


-1.56 


e  -1.565 
.2 
c 
o 


E 

re 


-1.57 


-1.575 


-1.58 


1 

j 

1 

j 

\  j 
\  i 

\/ 

V 

- sA 

V 

j 

50 


100 

Time  (days) 


150 


200 


Figure  10:  Hamiltonian  (Min  Time  to  Mars  Rendezvous) 


C.  MINIMUM  FUEL  RENDEZVOUS 

Next,  the  minimum  fuel  rendezvous  was  solved.  The  problem  was  set  up 
identically  to  the  minimum  time  problem  with  one  notable  exception,  the  initial  mass 
event  condition  was  free.  As  in  the  minimum  time  problem  with  initial  mass  free,  the 
solution  chose  to  depart  Earth  with  the  minimum  amount  of  mass  required  to  complete 
the  journey.  This  is  hardly  in  the  spirit  of  planetary  exploration,  therefore  the  cost 
function  was  modified  to  instead  maximize  the  final  mass  as  in  Eqn.  (96).  Thus  in  this 
case  “minimum  fuel”  is  a  bit  of  a  misnomer;  instead,  minimum  fuel  is  implied  by  the  true 
cost,  maximized  final  mass. 

A  second  difficulty  with  minimum  fuel  trajectories  is  that  they  tend  to  gravitate 
toward  extremely  long  trajectories  in  both  time  and  path  length.  If  both  time  and  transfer 
angle  are  unconstrained,  the  vehicle  will  begin  the  trajectory  with  a  very  short  duration 


45 


burn  which  slightly  raises  the  aphelion  of  the  trajectory.  The  vehicle  then  coasts  for 
slightly  more  than  a  year  before  returning  to  perihelion  where  it  again  conducts  a  short 
bum,  again  slightly  raising  aphelion.  This  pattern  is  repeated  over  many  years  until  the 
aphelion  has  finally  reached  the  target  planet’s  orbital  radius.  These  extremely  long 
trajectories  are  difficult  to  capture  using  a  direct  method  such  as  DIDO  due  to  the 
comparatively  few  number  of  nodes  typically  used.  To  prevent  this  from  occurring,  the 
transfer  angle  was  bound  such  that  only  type  I  orbits  would  be  permissible  (a  type  I  orbit 
is  an  orbit  with  a  total  transfer  angle  between  zero  and  180  degrees,  a  type  II  orbit 
between  180  and  360  degrees,  etc. )  That  is  to  say 

0  <  0  (f^ )  <  7C  (131) 

The  optimal  solution  presented  below  is  comprised  of  60  nodes  and 
required  7.85  minutes  to  converge.  The  total  flight  time  was  0.695694  years  (253.4  days) 
or  70  days  longer  than  the  minimum  time  trajectory.  However  only  119.1  kg  of 
propellant  are  used  giving  an  arrival  mass  of  540.3  kg.  Note  also  that  the  final  transfer 
angle  is  exactly  71  ,  the  maximum  allowable.  Once  again  the  propagated  states  (plotted  as 
lines)  are  in  excellent  agreement  with  the  DIDO  states  (plotted  as  various  symbols). 
Thus  the  trajectory  is  at  a  minimum  feasible. 
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3.5 


Low  Thrust  States:  Earth  to  Mars 


Figure  11:  State  History  (Min  Fuel  to  Mars  Rendezvous) 


The  heliocentric  view  of  the  trajectory  [Figure  12]  shows  how  this  mass  savings  is 
attained.  The  vehicle  begins  by  applying  a  near-tangential  burn  for  approximately  40 
days  before  shutting  off.  This  long  bum  places  the  vehicle  in  a  transfer  orbit  whose 
aphelion  exactly  corresponds  with  Mars’  orbit  radius.  The  vehicle  coasts  along  this 
trajectory  until  around  the  220th  day  where  it  begins  a  shorter,  sustained  burn  to  increase 
the  vehicle’s  speed  to  match  that  of  Mars. 
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Low  Thrust:  Earth  to  Mars 


Figure  12:  Heliocentric  Trajectory  (Min  Fuel  to  Mars  Rendezvous) 

The  next  plot  shows  the  control  history  throughout  the  trajectory.  Again  the  CMT 
and  DIDO  controls  match  extremely  well  for  the  thrust  angle  but  diverge  for  the  thrust 
magnitude  control.  The  DIDO  thrust  angle  control  (circles)  can  be  seen  to  become  more 
irregular  in  the  middle  portion  of  the  trajectory  when  vehicle  is  not  thrusting.  This  is 
expected  as  the  thrust  angle  is  not  well  defined  without  a  thrust  magnitude.  Much  more 
interestingly  is  the  fact  that  the  CMT  thrust  angle  control  history  (dots)  does  appear  to 
follow  a  well-defined,  smooth  curve  even  when  the  thrust  magnitude  is  zero.  Also  note 
that  although  the  thrust  magnitude  CMT  and  DIDO  controls  do  not  agree  during  the 
thrusting  portion  of  the  trajectory,  they  do  agree  during  the  non-thrusting  portion.  Thus  it 
seems  that  the  “switch”  is  only  working  in  one  direction. 
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Control  Histories 
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Let  us  investigate  this  discrepancy  further.  Recall  that  the  necessary  condition  for 
the  thrust  magnitude  is 


mag 


2(i 


rel 


X 

— -  +  Br 

y  mag 

V  e  J 


which  can  be  rearranged  as 


^Tmag^rel  + 


X. 


(132) 


(133) 


Because  DIDO  returns  both  the  covectors  with  the  solution  (via  the  CMT)  we  can  plug 
the  DIDO  solution  states,  controls  and  covectors  into  the  right  hand  side  of  Eqn.  (133) 
and  check  that  the  result  matches  the  thrust  magnitude  covector  as  returned  by  DIDO.  As 
we  can  see  in  the  following  plot,  the  two  do  not  agree. 
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Figure  14:  Comparison  of  DIDO  and  CMT  Thrust  Covector 
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We  are  left  with  two  possibilities,  either  Eqn.(133)  is  incorrect,  or  the  DIDO 
supplied  thrust  covector  is  incorrect. 

Despite  the  fact  that  we  can  not  verify  the  DIDO  controls  by  comparison 
to  the  CMT,  we  can  still  verify  the  thrust  magnitude  switches  by  plotting  the  switching 
function.  Application  of  the  KKT  theorem  tells  us  that 
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>0  if 
=  0 
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mag  max 

0  <  T  <T 

mag  m 


(134) 


This  implies  that  the  controls  switch  between  bounds  at  the  zero  crossings  of  the  covector 
history.  Figure  15  shows  the  thrust  magnitude  profile  plotted  above  the  switching 
function.  The  lower  plot  has  been  cropped  to  just  either  side  of  the  x  axis  to  better  see  the 
zero  crossings.  It  is  clear  that  the  zero  crossings  of  the  lower  plot  coincide  with  the  thrust 
switches  as  shown  in  the  upper  plot. 
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Thrust  and  Switching  Functions 
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Figure  15:  Switching  Function  (Min  Fuel  to  Mars  Rendezvous) 
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Finally,  the  Flamiltonian  is  plotted  in  Figure  16.  Once  again  the  Hamiltonian  is 
approximately  constant,  verifying  the  necessary  condition  (via  first  integral).  Note  that  in 
this  case  the  value  of  the  Hamiltonian  is  constant  at  zero.  This  is  expected  when  time  is 
not  explicitly  present  in  the  cost  function. 
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Figure  16:  Hamiltonian  (Min  Fuel  to  Mars  Rendezvous) 
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VII.  OPTIMAL  AEROCAPTURE  PROPLEM  FORMULATION 


A.  EVENT  CONDITIONS 

A  spacecraft  enters  a  planet’s  SOI  with  a  known  velocity  relative  to  the  planet. 
This  velocity,  known  as  arrival  V-infinity  ( )  determine  the  shape  of  the  trajectory  in 

the  planet  relative  frame.  When  >0  ,  as  in  the  case  of  aerocapture,  this  trajectory 

will  be  a  hyperbola  with  the  planet  at  the  focus.  For  trajectories  that  do  not  intersect  the 
atmosphere,  the  vehicle’s  trajectory  is  Keplerian  and  easily  solvable  at  any  point.  For  an 
aerocapture  trajectory,  the  vehicle  will  follow  such  a  Keplerian  trajectory  from  the  arrival 
point  to  the  upper  limit  of  the  atmosphere,  sometimes  referred  to  as  the  atmospheric 
interface.  The  atmospheric  limit  is  defined  as  the  altitude  above  the  planet’s  surface 
below  which  non-conservative  forces  such  as  atmospheric  drag  begin  to  perturb  the 
trajectory.  Thus  only  the  atmospheric  portion  of  the  aerocapture  trajectory  need  be 
optimized  for  two  reasons: 

■  All  portions  of  the  trajectory  outside  the  atmosphere  have  closed-form 
solutions  that  can  be  determined  uniquely  from  the  vehicle’s  velocity 
vector  evaluated  at  the  atmospheric  entry  point  (working  backward  to  V 
infinity  of  arrival)  and  atmospheric  exit  point  (working  forward  to  the 
apoapsis). 

■  Above  the  atmospheric  limit,  the  controls  have  no  affect  on  the  trajectory 
as  there  are  no  aerodynamic  forces 
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Atmosphere  Limit 
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Figure  17:  Atmospheric  Entry  and  Exit  Points 


This  being  the  case  we  can  formulate  our  initial  and  final  radius  event  conditions 
as 


(  ^0  )  ^atm  limit 

(135) 

(t  )  =  r 

\  f  )  atm  limit 

(136) 

As  was  the  case  with  initial  angular  displacement  in  the  low  thrust  problem,  the 
initial  longitude  is  invariant.  Thus  it  is  convenient  to  define  the  initial  longitude  (as 
defined  by  atmospheric  entry)  to  be  zero. 

e(f„)=o  (137) 
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The  well  know  vis -viva  equation  [Ref.  2]  defines  the  specific  energy  of  a  body  at 
any  point  outside  the  atmosphere  as  a  function  of  the  vehicle’s  inertial  speed  and  radial 
distance  from  the  central  body. 


2  r 


However  the  arrival  point  is  defined  as  r  — >  °o  so 

vl 


’'arrival 


Due  to  conservation  of  energy, 


£  =  £ 

arrival  atm— in 


so  substituting  Eqn.  (138)  and  (139)  into  Eqn.  (140)  we  get 


vt  y2  . 
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atm 


where  ratm  is  the  atmospheric  limit.  This  can  be  rearranged  to 


y2  -y2  .  I, 
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(138) 


(139) 


(140) 


(141) 


(142) 


Eqn.  (142)  becomes  the  equality  constraint  relating  the  vehicle’s  inertial  velocity 
at  atmospheric  entry  to  V ^  .  From  the  calculations  in  Appendix  A,  we  can  further 

relate  the  inertial  velocity  at  atmospheric  entry  to  the  states  as  measured  in  the  rotating 
frame  using 

y2  =v02  +  rQ2£l2  cos2  (f>0 +  2r0v0£2cos(|)0  cost)/  Ocoyy0  (143) 


where  all  state  variables  on  the  right-hand  side  are  taken  at  atmospheric  entry  (initial 
time,  also  recall  that  lower-case  velocity  state  symbols  represent  non  -inertial 
components). 

Similarly,  we  require  an  event  condition  to  ensure  that  the  vehicle’s  state  at 
atmospheric  exit  is  sufficient  to  target  our  desired  apoapsis  as  show  in  Figure  18  below. 
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Target  Circular 
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Figure  18:  Targeting  a  Circular  Orbit 


Again  due  to  conservation  of  energy  we  have 


£  =  £ 

atm— out  apo 


(144) 


and  using  vis -viva  we  obtain 


-17  2  Tf  2 

V  atm-e  u  t  |X  '  apo 

9  r  9  r 

^  'am  ^  'apo 


(145) 


so  V  is  needed  as  a  function  of  the  state  vector  at  atmospheric  exit. 


The  angular  momentum  of  a  point  mass  is  given  by 

h  =r  xV  =  rFcos(r)  (146) 

where  T  is  the  inertial  flight  path  angle.  Due  to  conservation  of  angular  momentum 
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(147) 


giving  us  an  equality  constraint  relating  the  states  at  atmospheric  exit  to  our  desired 
apoapsis.  Again,  using  the  relations  derived  in  Appendix  A  we  can  find  Vatm_out  and 

r utm  ou,  as  functions  of  the  state  variables  in  the  rotating  frame  using  the  following 
relations: 

V2  =  vf2  +r2Ct  cos2^  +  2rfvfQ.cosfyf  cost)/ /covy f  (152) 

v/sinY/ 

cos2y/  (v^.2  +f22r/2)  +  2v/r/f2cosy/  cost \rf  cost]^ 

where  in  this  case  all  the  state  variables  on  the  right-hand  side  are  again  taken  at 
atmospheric  exit. 

Note  that  these  event  conditions  require  that  the  apoapsis  of  the  post-atmospheric 
trajectory  exactly  intersect  the  circular  target  orbit,  also  known  as  a  Hohmann  transfer. 
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Hohmann  showed  that  a  Hohmann  transfer  is  the  minimum-fuel  transfer  orbit  between 
co-planar  orbits  whose  ratio  of  radii  are  less  than  1 1.94  [Ref.  7].  [Ref.  24]  develops  the 
event  conditions  necessary  for  a  bi-elliptic  transfer  orbit  where  the  apoapsis  of  the  post- 
atmospheric  trajectory  has  insufficient  energy  to  reach  the  circular  target  orbit.  Using  the 
same  approach,  one  could  likewise  develop  the  event  conditions  for  the  case  where  the 
post-atmospheric  trajectory  has  excess  energy  and  overshoots  the  target  circular  orbit. 
Unfortunately,  these  three  cases  (apoapsis  undershoot,  touching,  and  overshoot)  can  not 
easily  be  reconciled  into  one  set  of  conditions  because  doing  so  would  complicate  the  AV 
calculations  in  the  next  section,  necessitating  the  introduction  of  an  absolute  value 
operator  in  the  rocket  equation  Eqn.  (156).  Because  the  derivative  of  the  absolute  value 
operator  is  discontinuous  when  evaluated  at  zero,  singularities  are  introduced  in  the 
Jacobian  leading  to  serious  numerical  issues  in  solving  the  NLP.  [Ref.  5].  For  these 
reasons,  only  the  “touching”  case  is  considered  in  this  work. 

The  event  conditions  can  be  summarized  in  vector  form  as  as 
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r(ff) 
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0 
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(154) 


Where  Va,m-ou,  and  Tatm-ou,  given  bY  E9nS-  (152)  and  (153)' 


B.  COST  FUNCTIONS 

1.  Minimum  Fuel  to  Circularize 

The  simplest  performance  index  for  aerocapture  is  to  minimize  the  fuel  required  to 
circularize  the  orbit.  Since  the  event  conditions  require  the  post -atmospheric  trajectory’s 
apoapsis  exactly  touch  the  target  circular  orbit,  the  delta- V  can  be  obtained  by  subtracting 


the  apoapsis  velocity  found  using  Eqn.  (149)  from  the  target  circular  orbit  velocity 
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v„. 


(155) 


such  that 


AV=V.  -V 

circ  apo 


With  the  delta- V  known,  the  rocket  equation 


mprop  =  mo 


{  _AV  A 

/, 


(156) 


(157) 


l  —  e 

v  J 

can  be  used  to  find  the  propellant  mass  required  to  circularize  the  orbit.  Thus  our 
performance  index  can  be  expressed  in  Mayer  form  as 


J  — m „ 


(158) 


2.  Minimum  Heat  Load 

Knowing  that  there  exists  some  relation  between  the  total  heat  load  and  the 
required  TPS  mass,  another  reasonable  performance  index  is  to  minimize  the  heat  load. 
As  one  recalls,  heat  load  is  the  integrated  heating  rate  over  the  trajectory  (Eqn.  (50))  ; 
thus  a  natural  choice  is  to  use  a  Lagrange  form  cost  function  such  as 

'/ 

/  -  Jd(  A  v)  dx  (159) 

h 

However  any  Lagrange  cost  can  be  re-written  as  a  Mayer  cost  through  the 
addition  of  a  state  in  the  equations  of  motion.  In  this  case,  we  could  add  heat  load  as  a 
state,  with  dynamics  given  by  Eqn.  (50). 


3.  Minimum  Aerocapture  Mass 

A  more  useful  performance  index  would  to  minimize  the  total  mass  associated 
with  aerocapture.  The  total  aerocapture  mass  can  be  broken  into  three  components,  heat 
shield  mass  (also  known  as  fore-shield  mass),  back-shell  mass,  and  the  propellant  mass 
required  to  circularize  the  orbit. 
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^ heatshield  ^%nrk  —shell  ^ 


prop 


(160) 


The  propellant  mass  to  circularize  can  be  found  using  Eqn.  (157).  The  mass  of 
the  back-shell  can  be  assumed  to  be  relatively  constant  over  a  large  range  of  heat.  Thus  a 
method  is  needed  for  determining  the  mass  of  the  heat  shield. 

A  first-order  approximation  for  heat  shield  mass  can  be  made  by  linearly  mapping 
heat  load  to  heat  shield  mass.  Thermal  analysis  of  the  heat  loads  expected  to  be 
encountered  by  Mars  Smart  Lander  require  an  approximate  heat  shield  mass  of  40  kg  for 
every  10,000  Joules  of  heat  load  [Ref.  25].  Using  this  single  data  point  and  assuming  a 
20%  margin,  a  crude  mapping  between  heat  load  and  heat  shield  mass  can  be  expressed 
as 

>f 

mhea, shield  =  k\ <1  ^  (161) 

1 o 

where  k  =  50k§  . 

10,000 J 

Since  the  back-shield  mass  is  assumed  to  be  constant,  it  can  be  neglected  from  the 
performance  index.  A  minimum  aerocapture  mass  performance  index  can  now  be  written 
in  Bolza  form  as 


‘f 

J  =  m prop  +  k\q  dx  (162) 

It  should  be  stressed  that  a  much  more  accurate  model  could  be  obtained  by 
mapping  the  heat  shield  mass  to  both  heat  load  and  peak  heating  rate.  However,  the 
addition  of  the  peak  heating  rate  term  would  require  the  cost  function  to  be  formulated  as 
a  Chebyshev  minimax  problem  which  is  beyond  the  scope  of  this  proof-of-concept. 
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c. 


CONTROLS 


For  non-thru sting,  fixed  angle-of-attack  aerocapture  trajectories,  the  only  control 
during  the  atmospheric  pass  is  the  bank  angle  of  the  spacecraft.  As  in  the  low  thrust  case, 
coding  the  control  as  an  angle  and  bounding  it  as  0  <  8  <  271  leads  to  some  difficulties 
due  to  the  jump  discontinuity  between  0  and  2k.  As  shown  in  later  sections,  the  optimal 
bank  angle  history  generally  switches  between  7t  and  2rc  (or  zero).  Opening  the  control 
bound  to  0  <8  <371  allows  for  optimal  bank  angles  that  are  entirely  interior.  However  as 
in  the  low  thrust  case  it  was  eventually  realized  that  using  sin8  and  cos8  as  controls 
tended  to  yield  faster  and  more  consistently  reliable  solutions.  Thus  we  have 


(163) 


where  our  controls  are  bound  by 


-1  <  sin  8  <1 

(164) 

-1  <cosS  <  1 

(165) 

As  in  the  low  thrust  case,  this  formulation  of  the  controls  requires  the  addition  of 
a  path  constraint  (given  in  the  next  section). 

An  alternate  control  scheme  incorporates  the  bank  angle  as  a  state  variable  of  the 
vehicle  and  instead  controls  the  bank  angle.  This  allows  for  more  realistic  command 
response  behavior.  However,  this  scheme  was  not  incorporated  into  this  work. 


D.  STATE  BOUNDS 

Next  suitable  choices  must  be  made  to  bound  the  states.  Since  only  the 
atmospheric  portion  of  the  trajectory  is  of  concern,  there  is  no  point  in  allowing  the  radius 
to  be  any  larger  than  the  atmospheric  limit.  Similarly,  the  vehicle  has  a  hard  lower 
bound  at  the  planet’s  surface.  Considering  this  yields 


<  r<r. 


(166) 


However,  recalling  Eqn.  (49),  the  atmospheric  model  is  only  valid  to  the  scale 
height  below  which  the  exponential  atmosphere  incorrectly  decreases  with  decreasing 
altitude.  Thus  a  better  set  of  bounds  is 
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(167) 


planet 


atmos-limit 


Realizing  that  for  a  non-thrusting  vehicle  entering  the  atmosphere,  the  maximum 
speed  will  be  near  that  at  atmospheric  interface.  To  account  for  the  unlikely  event  of 
extremely  steep  entries  the  upper  bound  can  be  set  to  some  arbitrary  small  multiple  of  the 
initial  velocity. 


0<  v<5v 


atmos  -in 


(168) 


Note  however  that  Eqns.  (25),  (26),  (31),  (33)  have  singularities  at  v  =  0.  This  can  be 
resoled  by  constraining  the  problem  as 

£  <v<v,  (169) 

v  —  —  atmos- m  / 


where  ev  is  number  slightly  larger  than  zero. 

Note  that  equations  Eqns.  (22),  (25),  (31)  also  have  singularities  at  <|)  =  ±Jt  and 
y  =±n;  .  We  can  prevent  the  singularities  from  being  encountered  by  bounding  these 
states  as 


-n+%  ^  ^  -e„ 

(170) 

—71  +£y<  y<7t  -£y 

(171) 

where  again  £0  and  £y  are  small  numbers  greater  than  zero. 


Finally  consider  longitude  and  heading  angle.  Since  there  are  no  singularities  or 
physical  bounds  to  contend  with  we  are  free  to  choose  any  convenient  set  of  bounds 
sufficiently  large  to  reduce  the  NLP  search  space.  Note  however  that  in  both  cases 


0< 


0 

¥ 


<  271 


(172) 


is  a  poor  choice.  While  we  as  humans  intuitively  know  that  these  variables  are 
continuous  across  these  bounds,  the  dynamics  are  blind  to  this  fact  and  instead  see  this  as 
a  discontinuity.  This  discontinuity  can  be  resolved  by  simply  opening  the  bounds 
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The  bounds  for  the  complete  state  vector  can  be  summarized  as 
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E.  PATH  CONSTRAINTS 

The  control  space  in  inz9l2  defined  by  Iqns.  (164)  and  (165)  is  the  area  of  a 
square.  A  path  constraint  is  needed  to  correct  this  control  space  to  that  of  a  circle.  This 
is  accomplished  using  the  trigonometric  identity 

sin2  8  +cos2S  -1  =  0  (175) 

However  this  control  surface  leads  to  a  convexity  problem.  Using  the  same 
technique  employed  for  the  low  thrust  control,  the  path  constraint  is  modified  to  an 
inequality  constraint 

-1  <  sin25 +cos25 -1  <  0  (176) 

This  transforms  the  control  space  to  a  circular  disc  of  unit  radius  which  is  a 
convex  surface.  Again,  this  temporary  violation  of  physics  will  be  resolved  by  the 
Minimum  Principle  which  forces  the  controls  to  the  boundary  of  the  control  space,  in  this 
case,  onto  the  unit  circle. 

Other  path  constraints  of  potential  importance  to  the  mission  designer  are  limiting 
maximum  dynamic  pressure,  load-factor  or  heating-rate  as  given  by  Eqns.  (37),  (42)  and 
(50)  respectively. 
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F. 


NECESSARY  CONDITIONS 

Assuming  a  Mayer  cost  (Fe  0 ),  the  Hamiltonian  for  the  aerocapture  problem  is 


H  =  Xr  (vsiny)  +  A,e 


f  v cosy  cost)/  ^ 
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Note  that  all  of  the  aerocapture  Lagrange  costs  are  pure  state  costs  meaning  that  they  are 
not  function  of  the  controls.  Since  the  partial  of  the  augmented  Hamiltonian  will  be  taken 
in  a  moment,  the  assumption  of  a  Mayer  cost  does  not  impact  the  formulation  of  the 
necessary  conditions. 

The  augmented  Hamiltonian  is 


H'=\ 


LsinS 
mv cosy 


+  \ 


Lc  os  8 


v  mv  j 


+  M*  sin5 


(178) 


+ 


|lc5  cosS  +(ire/  (sin28  +cos2S  -l)  +  pure  state  terms 


where  all  the  pure-state  terms  have  been  grouped  together  and  where  (ls8 ,  (lf5  and  jlrel 
are  the  Lagrange  multipliers  associated  with  the  two  controls  and  the  path  constraint 
relating  them  respectively  [Eqn.  (176)]. 

Applying  the  necessary  condition  for  optimality  by  taking  the  partial  derivative  of 
the  augmented  Hamiltonian  with  respect  to  each  control  we  obtain 
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(179) 


dH'  _\L 
d  (cos  8)  mv 


+  (Ls  +2|ire,cos8  =  0 


(180) 


which  can  be  re-arranged  to  solve  for  the  controls 
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sinS  = - ^ - (181) 

2m  cosy  2\irel 

cos  8  = - — - (182) 

2mv[lrel  2jirel 
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Yin.  OPTIMAL  AEROCAPTURE  RESULTS 


A.  PROCESS 

Problems  were  solved  beginning  with  easier  problem  formulations  before  moving 
on  to  more  challenging  cost  functions  and  higher  fidelity.  Solutions  were  obtained  for 
non-rotating  planets  with  fixed  initial  conditions  and  minimum  fuel  to  circularize  as  the 
performance  index.  Next  minimum  heat -load  trajectories  were  solved  before  attempting 
rotating  atmospheres.  With  these  successes  in  hand,  the  initial  conditions  were  relaxed 
such  that  the  velocity  at  atmospheric  entry  need  only  agree  with  the  arrival  V- infinity. 
This  allowed  for  solutions  with  excess  Vinfinity  at  arrival  to  be  obtained.  Finally,  the 
solution  for  minimum  total  aerocapture  mass  was  obtained. 


B.  MINIMUM  AEROCAPTURE  MASS  AT  MARS  WITH  ZERO  ARRIVAL 

V-INFINITY 

This  case  considered  a  vehicle  arriving  at  a  rotating  Mars  with  zero  excess 
velocity  with  minimum  total  aerocapture  mass  as  the  performance  index.  This  initial 
arrival  condition  could  occur  in  the  case  of  a  low  thrust  heliocentric  trajectory  whereby 
the  interplanetary  trajectory  is  a  rendezvous. 

Although  only  30  nodes  were  necessary  to  obtain  a  solution,  using  50  nodes 
greatly  increased  the  accuracy  of  the  solution.  The  resulting  performance  index 
breakdown  is  given  in  Table  4.  The  total  aerocapture  mass  (minus  back-shell  mass  which 
is  assumed  constant)  is  only  19.14  kg.  This  low  value  can  be  attributed  to  poor  modeling 
(poor  choice  for  k )  in  Eqn.  (161)  as  well  as  the  low  heat  loads  experienced  due  to  the 
zero  V-infinity  at  arrival. 


Propellant  mass 

10.63  kg 

Front-shield  mass 

8.5  kg 

Total 

19.14  kg 

Table  4:  Cost  Function  Breakdown  (Zero  Arrival  V-inf) 
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Even  taking  these  factors  into  account,  the  optimal  aerocapture  mass  is  of 
significant  savings  when  compared  to  the  required  mass  for  a  pure  propulsive  injection 
maneuver.  Using  the  same  main  engine  and  bi-propellant  modeled  for  the  post- 
aerocapture  circularization  maneuver,  the  required  propellant  mass  for  the  pure- 
propellant  capture  is  201.44  kg.  Figure  19  shows  the  altitude,  latitude  and  longitude 
histories  during  the  atmospheric  portion  of  the  trajectory.  The  circles  represent  the  DIDO 
solutions  at  the  node  points,  whereas  the  line  represents  the  Runga-Kutta  propagated 
solution  (note  the  strong  agreement  between  the  two  indicating  feasibility).  As 
constrained,  the  trajectory  begins  and  ends  at  the  defined  atmospheric  interface  of  125 
km.  The  minimum  altitude  of  70.14  km  occurs  at  the  2.65  minute  mark  of  the  12.24 
minute  trajectory.  Note  that  the  initial  latitude  solution  for  the  pass  is  very  nearly  zero 
(0.7  deg)  and  that  the  trajectory  proceeds  easterly. 
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Figure  19:  Position  History  (Zero  Arrival  V-inf) 


Figure  20  shows  the  result  of  propagating  the  trajectory  beyond  the  atmospheric 

interface.  The  spacecraft  proceeds  to  an  altitude  of  295.4  km,  only  4.7  km  in  error  of  the 
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targeted  altitude.  This  error  is  attributed  to  the  error  in  the  final  states  at  atmospheric  and 
will  be  discussed  more  fully  ahead. 


Propagated  to  Apoapsis 


Figure  20:  Propagated  to  Apoapsis  (Zero  Arrival  V-inf) 


The  velocity  state  histories  are  given  in  Figure  21  and  demonstrate  the 
“backwards-S  curve”  in  speed  characteristic  to  aerocapture  trajectories.  The  trajectory 
begins  with  an  inertial  speed  of  4949.3  m/s  and  a  relatively  shallow  flight  path  angle  of 
-7.4  deg.  The  initial  heading  is  only  -0.09  deg,  agreeing  with  the  due  easterly  track  in. 
Taken  together,  it  is  clear  that  the  optimal  solution  is  to  fly  in  the  direction  of  the  rotating 
atmosphere  at  location  of  the  atmosphere’s  greatest  velocity  (the  equator).  This 
trajectory  gives  the  minimum  relative  speed  between  the  atmosphere  and  the  vehicle,  thus 
reducing  the  heating-rate  which  is  proportional  to  VM.  It  is  interesting  to  note  that  even 
when  the  problem  is  seeded  with  a  westerly-tracking  guess,  DIDO  still  returns  an  easterly 
optimal  solution.  The  inertial  velocity  at  atmospheric  exit  is  3523.5  m/s  for  a  total 
aerocapture  delta- V  of  1425.8  m/s. 
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Velocity  History 
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Figure  21:  Velocity  History  (Zero  Arrival  V-inf) 


The  accuracy  of  the  solution  can  be  evaluated  by  comparing  the  DIDO  terminal 
states  with  those  from  the  propagator  as  shown  in  Table  5.  The  fact  that  the  altitude 
indeed  hits  its  target  value  (the  only  explicitly  constrained  final  state)  demonstrates  the 
accuracy  of  the  solution.  The  error  that  remains  is  likely  due  to  these  errors  being  less 
than  the  tolerances  set  in  DIDO,  thus  causing  the  solver  to  exit.  Furthermore,  radius,  not 
altitude  is  the  actual  state  vector  used  in  the  solution.  Redefining  the  equations  of  motion 
to  use  altitude  as  a  state  instead  may  allow  for  increased  accuracy. 
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State 

DIDO  value 

Propagator  value 

Absolute  Error 

Percent  Error 

Hf 

125  km 

124.52  km 

-0.4762  km 

0.3824  % 

e/ 

44.68  deg 

44.67  deg 

-0.0111  deg 

0.0248  % 

4 >/ 

0.11  deg 

0.11  deg 

0.0005  deg 

0.4335  % 

vf 

3274.5  m/s 

3272.7  m/s 

-1.8559  m/s 

0.0567  % 

¥/ 

0.09  deg 

0.09  deg 

0.0011  deg 

1.2560  % 

Y/ 

2.02  deg 

1.98  deg 

-0.0328  deg 

1.6523  % 

Table  5:  Propagated  Accuracy  (Zero  Arrival  V-inf) 


Figure  22  shows  the  optimal  bank  angle  history  during  the  aerocapture  pass.  The 
blue  circles  represent  the  DIDO  solution  controls,  whereas  the  red  dots  represent  the 
CMT  derived  controls  evaluated  using  Eqns.  (181)  and  (182).  For  the  bulk  of  the 
trajectory,  the  bank  angle  takes  one  of  two  approximate  values,  0  deg  (lift-up)  or  180  deg 
(lift -down).  This  choice  of  extreme  controls  often  occurs  in  optimal  control  solutions.  At 
some  points,  the  DIDO  control  solution  actually  oscillates  between  -180  deg  and  180  deg 
which  of  course  are  equivalent.  The  same  controls  have  been  shifted  to  lie  between 
0<5  <271  in  Figure  23  to  reduce  this  numerical  irregularity.  Again,  note  the  excellent 
agreement  between  the  DIDO  solution  controls  and  the  CMT  controls.  The  only 
significant  divergence  between  the  two  occurs  late  in  the  problem  when  the  vehicle  is 
near  its  slowest  and  highest  portion  of  the  trajectory.  In  this  energy  state,  the 
performance  index  is  of  greatly  reduced  sensitivity  to  the  control  (in  other  words,  these 
errors  are  unimportant  to  the  solution.)  This  agreement  verifies  that  the  first-order 
necessary  conditions  have  been  met  strengthening  the  argument  that  the  solution  is  at 
least  a  local  minimum.  Furthermore,  the  flatness  of  the  Hamiltonian  in  Figure  24  proves 
verification  of  the  first  integral  further  contributing  to  the  claim  of  optimality. 
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Figure  22:  DIDO  and  CMT  Controls  (Zero  Arrival  V-inf) 
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Figure  23:  DIDO  and  CMT  Controls  (0-2  % )  (Zero  Arrival  V-inf) 
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The  bank  angle  “switch”  occurs  at  approximately  3.2  minutes,  fully  30  seconds 
after  minimum  altitude  passage.  A  possible  explanation  for  this  switching  strategy  is  as 
follows:  a  steep,  lift-up  entry  trajectory  allows  for  deep  atmospheric  penetration  with 
positive  curvature.  Once  the  “corner  has  been  turned”  at  nadir  and  the  velocity  state  is 
such  that  atmospheric  exit  is  guaranteed,  the  vehicle  flips  lift-down  essentially  “clinging” 
to  the  atmosphere.  This  extends  the  time  within  the  atmosphere  but  on  the  slower  side  of 
the  minimum  altitude  where  the  heating  rate  penalty  is  less  severe. 
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Figure  24:  Hamiltonian  (Zero  Arrival  V-inf) 


The  heating  rate  over  the  trajectory  is  given  in  Figure  25.  The  peak  heating  rate 
of  7.44  W/cm  occurs  at  140.18  sec  (2.34  min).  The  total  heat  load  was  calculated  as 

'y 

1700.38  J/cmr  from  the  heating  rate  history  using  a  trapezoidal  integration  scheme 
(MATLAB’s  “ trapz ”  command). 
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Stagnation  Point  Heating 
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Figure  25:  Heating  rate  (Zero  Arrival  V-inf) 


The  following  figures  show  the  time  history  of  the  body  accelerations  resolved  in 
the  Frenet  frame  and  shows  the  total  acceleration  peaking  at  -0.85  g’s  at  time  158.94  sec 
(2.65  min),  the  vast  majority  of  which  is  in  the  anti-tangential  direction  with  small  normal 
and  bi-normal  components.  Not  surprisingly,  the  peak  dynamic  pressure  (shown  in 
Figure  28)  of  505.71  Pa  occurs  at  the  same  time  (recall  from  Eqn.  (38)  that  the  tangential 
acceleration  is  a  function  of  drag  and  mass  only). 
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Figure  26:  Body  Accelerations  (Zero  Arrival  V-inf) 
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Figure  27:  Total  Acceleration  (Zero  Arrival  V-inf) 
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Dynamic  Pressure  vs  Time 


Figure  28:  Dynamic  Pressure  (Zero  Arrival  V-inf) 


Finally,  selected  orbit  elements  are  plotted  versus  time  in  Figure  29.  One  can  see 
that  the  non-dimensional  specific  energy  drops  from  zero  energy  (recall  that  the  vehicle 
begins  with  zero  arrival  V-infinity)  and  that  the  eccentricity  begins  at  a  parabolic  value  of 
one  before  decaying  to  a  near-circular  value  of  0.04  at  atmospheric  exit.  The  apoapsis 
begins  at  positive  infinity  (as  it  is  singular  for  a  parabola)  and  decays  to  a  non- 
dimensional  value  corresponding  to  300.0  km. 


76 


Orbit  Paramters  (Non-Dimensional) 


Time  (min) 

Figure  29:  Selected  Orbit  Parameters  (Zero  Arrival  V-inf) 


C.  MINIMUM  AEROCAPTURE  MASS  AT  MARS  WITH  EXCESS 

ARRIVAL  V-INFINITY 

This  case  considered  a  vehicle  arriving  at  a  rotating  Mars  with  an  excess  velocity 
with  minimum  total  aerocapture  mass  as  the  performance  index.  In  this  case,  the  initial 
arrival  velocity  is  5.2  km/s.  This  roughly  equates  to  the  arrival  velocity  expected  from  an 
impulsive  Hohman  transfer  between  Earth  and  Mars. 

The  solution  given  below  was  obtained  using  the  solution  from  Section  B  above 
and  bootstrapping  up  by  increments  of  10  nodes  until  a  satisfactory  solution  was  obtained 
with  90  nodes.  The  cost  breakdown  for  this  solution  is  given  in  Table  6.  Although  this 
trajectory  only  requires  one  kilogram  more  of  propellant,  the  required  heat  shield  mass  is 
slightly  more  than  double  the  amount  required  for  the  zero  arrival  V- infinity  solution.  As 
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shall  be  seen,  this  is  due  to  the  significantly  higher  thermal  loads  placed  on  the  vehicle 
during  the  atmospheric  pass. 


Propellant  mass 

11.81  kg 

Front-shield  mass 

16.94  kg 

Total 

28.75  kg 

Table  6:  Cost  Function  Breakdown  (Excess  Arrival  V-inf) 


For  this  case,  the  cost  savings  when  compared  to  a  pure-propulsive  capture  are 
even  more  significant  with  the  pure-propulsive  capture  requiring  a  whopping  386.4  kg. 
Figure  30  shows  the  altitude,  latitude  and  longitude  histories  during  the  atmospheric 
portion  of  the  trajectory.  Note  once  again  the  strong  agreement  between  the  DIDO  state 
history  and  the  propagated  solution.  The  minimum  altitude  of  58.5  km  occurs  94.63  sec 
(1.58  min)  into  the  601.02  sec  (10.02  min)  trajectory.  Again  the  optimal  solution  begins 
near  zero  latitude  and  tracks  easterly  with  the  rotating  atmosphere. 
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Position  History 


Figure  30:  Position  History  (Excess  Arrival  V-inf) 


Figure  31  shows  the  continued  propagation  beyond  the  atmospheric  exit  point.  In 
this  case  the  propagated  apoapsis  reaches  only  277.28  km,  22.72  km  short  of  the  targeted 
300  km  altitude.  This  error  is  likely  due  altitude  error  at  atmospheric  interface  (discussed 
shortly)  magnified  with  time. 
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Propagated  to  Apoapsis 


Figure  31:  Propagating  to  Apoapsis  (Excess  Arrival  V-inf) 


Figure  32  gives  the  velocity  state  histories  with  an  inertial  entry  velocity  of 
7178.8  km/s  and  initial  inertial  flight  path  angle  of  -10.47  deg.  The  exit  velocity  is 
3517.1  km/s  for  a  significant  delta- V  of  3661.8  m/s. 
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Velocity  History 
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Figure  32:  Velocity  History  (Excess  Arrival  V-inf) 
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As  shown  in  Figure  33,  this  solution’s  altitude  trajectory  begins  with  a  steeper 
flight  path  angle  and  penetrates  deeper  into  the  atmosphere  than  the  slower,  zero  arrival 
V-infinity  solution.  Note  that  although  it  penetrates  deeper,  the  time  of  passage  is 
significantly  shorter.  As  we  shall  soon  see,  this  contributes  to  a  larger,  but  narrower 
heating-rate  pulse  which  helps  to  reduce  the  total  heat  load. 
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Altitude  vs.  Time  Comparison 


Figure  33:  Altitude  Profile  Comparison  (Excess  Arrival  V-inf) 


The  accuracy  and  feasibility  of  the  solution  is  shown  in  the  following  table.  The 
error  between  the  propagated  solution  and  the  DIDO  solution  is  small  enough  to  declare 
convergence.  Note  that  although  the  final  flight  path  angle  has  a  large  percentage  error, 
the  absolute  error  is  less  that  0.1  deg.  Because  the  final  value  of  the  flight  path  is  itself 
small,  this  makes  the  error  appear  large. 
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State 

DIDO  value 

Propagator  value 

Absolute  Error 

Percent  Error 

Hf 

125  km 

123.67  km 

-1.3323  km 

1.0773  % 

e/ 

39.42  deg 

39.39  deg 

-0.0301  deg 

0.0763  % 

4 >/ 

0.20  deg 

0.20  deg 

0.0011  deg 

0.0763  % 

vf 

3268.11  m/s 

3262.32  m/s 

-5.7845  m/s 

0.1773  % 

¥/ 

-0.53  deg 

-0.53  deg 

0.0001  deg 

0.0270  % 

Y/ 

2.35  deg 

2.25  deg 

-0.0989  deg 

4.3985  % 

Table  7:  Propagated  Accuracy  (Excess  Arrival  V-inf) 


The  bank  control  history  is  given  in  Figure  34.  The  “bank-bang”  nature  of  the 
control  can  be  seen  as  the  vehicle  begins  the  trajectory  lift -up  before  abruptly  switching 
to  lift -down  at  approximately  1.8  minutes.  Again,  the  DIDO  solution  can  be  seen 
flipping  back  and  forth  between  +180  and  - 180  which  are  of  course  equivalent.  Both  the 
switch  as  well  as  the  agreement  between  the  DIDO  controls  and  the  CMT  controls  is 
more  prevalent  in  this  case,  likely  do  to  the  significantly  larger  number  of  nodes  in  the 
solution.  Once  again  this  satisfies  the  first  order  necessary  conditions  for  optimality. 
Again  the  switch  occurs  after  (but  closer  to)  minimum  altitude  passage  by  approximately 
13  seconds.  Hie  optimality  of  the  solution  is  further  supported  by  the  constancy  of  the 
first  integral  in  Figure  35  (that  is  to  say  the  flatness  of  the  Hamiltonian  to  the  10"  ). 
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Figure  34:  DIDO  and  CMT  Control  Histoiy  (Excess  Arrival  V-inf) 
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Figure  35:  Hamiltonian  (Excess  Arrival  V-inf) 
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The  heating  rate  at  the  stagnation  point  is  given  in  the  following  figure.  The  peak 
heating  rate  of  35.02  W/cnr  occurs  at  time  72.76  sec  (1.21  min)  and  the  total  heat  load  is 
3388.37  J/cm  -  almost  double  the  value  for  the  earlier  non-rotating  case.  As  mentioned 
previously,  the  duration  of  the  heat  pulse  is  quite  short,  less  than  a  minute  as  measured  at 
the  half-maximum  point.  This  contrasts  with  the  non-rotating  case  where  the  half¬ 
maximum  pulse  width  was  approximately  three  minutes  in  duration. 
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Stagnation  Point  Heating 
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Figure  36:  Stagnation  Point  Heating  Rate  (Excess  Arrival  V-inf) 


The  body  accelerations  are  significantly  stronger  than  those  experienced  in 
the  previous  case  [Ref.  Figure  37].  The  tangential  acceleration  peaks  at  -4.5  g’s 
while  the  normal  accelerations  peak  at  0.82  g’s.  The  total  acceleration  peaks  at 
4.57  g’s  at  time  79.78  sec  (1.33  min)  [Figure  38].  Again,  this  corresponds  exactly 
to  the  dynamic  pressure  peak  of  2703.70  Pa  [Figure  39]. 
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Body  Accelerations 


Acceleration  Mag  vs  Time 
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Figure  38:  Total  Acceleration  (Excess  Arrival  V-inf) 
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Dynamic  Pressure  vs  Time 


Figure  39:  Dynamic  Pressure  (Excess  Arrival  V-inf) 

The  orbit  parameters  in  the  following  figure  are  slightly  more  interesting 
than  in  the  previous  case.  With  this  solution,  the  non-dimensional  energy  can  be 
seen  to  begin  well  above  zero,  with  capture  occurring  at  110.50  sec  (1.84  min). 
The  eccentricity  decays  from  a  hyperbolic  3.16  to  a  near-circular  0.04.  Also  the 
fact  that  apoapsis  is  undefined  for  parabolic  trajectories  is  demonstrated  as  the 
apoapsis  departs  toward  negative  infinity  to  the  left  of  the  singularity  before 
rapidly  falling  off  from  positive  infinity  immediately  following  the  capture. 
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Orbit  Paramters  (Non-Dimensional) 


Figure  40:  Selected  Orbit  Parameters  (Excess  Arrival  V-inf) 


D.  MINIMUM  AEROCAPTURE  MASS  AT  NEPTUNE  WITH  EXCESS 

ARRIVAU  V-INFINITY 

The  next  case  considered  a  minimum  total  aerocapture  mass  at  Neptune  with 
excess  arrival  V-infinity.  The  arrival  V-infinity  at  Neptune  was  9.42  km/s,  significantly 
higher  than  that  of  Mars.  A  higher  altitude  of  1000  km  was  targeted  for  the  final  circular 
orbit  and  the  atmospheric  interface  was  defined  as  an  altitude  of  800  km. 

The  problem  was  solved  by  first  solving  the  zero  excess  arrival  trajectory  for  a 
non-rotating  atmosphere.  This  solution  was  used  to  bootstrap  the  excess  arrival  velocity 
which  was  in-turn  used  to  bootstrap  the  case  of  a  rotating  atmosphere.  90  nodes  were 
deemed  sufficient  to  provide  an  accurate,  optimal  solution.  Table  8  gives  the  cost 
function  breakdown  for  this  solution.  The  total  required  aerocapture  mass  is  328  kg  of 
which  48.4  kg  is  propellant  mass  and  279.5  kg  is  heat  shield  mass.  While  significantly 
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higher  than  the  aerocapture  mass  required  at  Mars,  the  total  capture  mass  is  200  kg  less 
than  the  528.7  kg  of  propellant  required  for  a  purely  propulsive  capture  maneuver. 
Again,  the  accuracy  of  the  heat  shield  mass  predictions  is  limited  to  the  model  used, 
which  assumes  heat  shield  mass  scales  linearly  with  heat  load. 


Propellant  mass 

48.4  kg 

Front-shield  mass 

279.5  kg 

Total 

328.0  kg 

Table  8:  Cost  Function  Breakdown  (Neptune  Excess  V-inf) 


Figure  41  shows  the  position  history  of  the  spacecraft  during  the  atmospheric 
portion  of  the  trajectory.  Note  the  excellent  agreement  between  DIDO  solution  (circles) 
and  the  propagated  solution  (solid  line).  The  trajectory  begins  and  ends  at  800  km  of 
altitude  (the  defined  atmospheric  interface).  The  total  pass  requires  22.64  minutes  with  a 
minimum  altitude  of  256.5  km  occurring  4.16  minutes  into  the  trajectory.  Again  the 
solver  chooses  an  initial  latitude  of  nearly  0  degrees  (equatorial)  with  an  easterly  track. 
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Position  History 


5  400h 
200j- 


0  DIDO 

—  Propagator 

— 

S’-O.l 
3 -0.2 
i°-3< 

-0.4 


U) 

"§ 

a> 

~o 

3 

'5> 

c 

o 


10 


15 


20 


25 


- 

10  15 

Time  (min) 


20 


25 


^eeeeeee^ 

■mbSSSS*®* 

10  15 

Time  (min) 


20 


25 


Figure  41:  Position  History  (Neptune  Excess  V-inf) 

Figure  42  shows  the  result  of  propagating  the  DIDO  solution  beyond  the 
atmospheric  limit  of  the  solution.  The  propagated  apoapsis  of  976.12  km  is  within  23.9 
km  of  the  targeted  1000  km  apoapsis  altitude. 
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Propagated  to  Apoapsis 


Figure  42:  Propagated  to  Apoapsis  (Neptune  Excess  V-inf) 


Figure  43  gives  the  velocity  state  histories  for  the  Neptune  solution.  The 
trajectory  begins  at  a  staggering  inertial  atmospheric  entry  speed  of  nearly  25000  m/s,  a 
flight  path  angle  of  -9.35  degrees,  and  heading  -0.1  degrees  (due  east).  Inertial 
atmospheric  exit  occurs  at  16144  m/s  for  an  aerocapture  delta-V  of  8844  m/s. 
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Figure  43:  Velocity  History  (Neptune  Excess  V-inf) 

The  accuracy  of  the  solution  is  presented  in  Table  9  which  compares  the  DIDO 
and  propagated  terminal  state  values.  Note  that  the  “large”  percentage  errors  in  latitude 
and  flight  path  angle  are  due  to  the  small  absolute  value  of  the  state. 


State 

DIDO  value 

Propagator  value 

Absolute  Error 

Percent  Error 

Hf 

800  km 

793.07  km 

-6.9257  km 

0.8733  % 

9/ 

48.32  deg 

48.31  deg 

-0.0116  deg 

0.0241  % 

<1 >/ 

-0.00  deg 

-0.00  deg 

-0.0007  deg 

19.5726  % 

vf 

13736.80  m/s 

13732.23  m/s 

-4.5689  m/s 

0.0333  % 

V/ 

0.64  deg 

0.64  deg 

-0.0007  deg 

0.1067  % 

Y/ 

1.51  deg 

1.46  deg 

-0.0495  deg 

3.3939  % 

Table  9:  Propagated  Accuracy  (Neptune  Excess  V-inl) 
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The  control  history  for  the  trajectory  is  presented  in  Figure  44.  The  DIDO 
solution  for  the  first  minute  and  last  nine  minutes  appears  somewhat  erratic. 
Investigating  this  further,  the  covectors  associated  with  the  controls  over  these  time 
intervals  are  of  very  small  value,  indicating  a  lack  of  sensitivity  of  the  performance  index 
to  the  bank  angle  in  these  regions.  This  corresponds  with  the  physical  explanation  of 
reduced  control  effectiveness  in  the  thinner  upper  limits  of  the  atmosphere.  In  the 
thicker,  lower  atmosphere,  the  bank  angle  trajectory  assumes  the  previously  seen,  lift- 
up/lift -down  “bang-bang”  type  control  with  the  switch  occurring  at  approximately  5.5 
minutes,  slightly  after  the  minimum  altitude  point.  The  DIDO  controls  (circles)  and  the 
CMT  derived  controls  (dots)  are  again  in  excellent  agreement. 
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Figure  44:  DIDO  and  CMT  Control  History  (Neptune  Excess  V-inf) 


The  thicker  atmosphere  and  significantly  atmospheric  velocities  contribute  to 
extreme  stagnation  point  heating  for  this  trajectory.  The  maximum  heating  rate  of 
264.59  W/cm  occurs  just  prior  to  the  minimum  altitude  point  (Figure  45).  The  long  time 
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duration  of  the  pass  results  in  a  total  heat  load  of  55,904  J/cm2  driving  up  the  mass 
requirements  for  the  TPS. 

Stagnation  Point  Heating 
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Figure  45:  Stagnation  Point  Heating  Rate  (Neptune  Excess  V-inf) 


Thermal  protection  is  not  the  only  engineering  challenge  presented  by  aerocapture 
at  Neptune.  The  large  aerodynamic  forces  lead  to  a  peak  total  acceleration  of  5.83  g 
(Figure  46).  These  high  loads  more  akin  to  those  experienced  by  fighter  aircraft  would 
require  additional  structural  mass,  further  reducing  available  payload  mass.  This 
additional  mass  cost  is  not  included  in  the  modeling  of  this  work. 
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Acceleration  Mag  vs  Time 


Figure  46:  Total  Acceleration  (Neptune  Excess  V-inf) 

E.  MINIMUM  AEROCAPTURE  MASS  AT  MARS  SUBJECT  TO  G-UIMITS 

The  previous  results  at  Neptune  demonstrate  the  potential  for  high  g-loads  during 
aerocapture  maneuvers.  Rather  than  increasing  the  structural  integrity  of  the  spacecraft 
to  survive  the  loads,  another  option  is  to  simply  constrain  the  g-limit  to  some  smaller, 
more  manageable  value.  The  following  solution  shows  partial  results  for  the  problem 
described  by  VIII. C  (Minimum  total  aerocapture  mass  at  Mars  with  excess  arrival 
velocity)  except  that  the  tangential  acceleration  has  been  limited  to  3  g.  A  120  node 
solution  was  obtained  by  bootstrapping  from  the  unconstrained  solution  results.  The 
trajectories  were  similar  with  the  constrained  trajectory  being  slightly  shallower  (initial 
flight  path  angle  of  -9.73  degrees  compared  with  -10.47)  and  longer  in  duration  (total 
pass  time  of  11.4  minutes  compared  with  10.02  minutes).  The  position  state  histories  are 
presented  in  Figure  47.  The  minimum  altitude  varies  by  only  3.5  km  but  occurs  almost 
30  seconds  later  than  the  unconstrained  case. 
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Position  History 
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Figure  47:  Position  History  (G-limited) 

The  velocity  histories  are  nearly  identical,  with  a  delta- V  within  4  m/s  of  the 
unconstrained  case.  The  only  notable  differences  are  some  slight  perturbations  in  the 
flight  path  angle  and  heading  angle  between  times  1.5  minutes  and  2.5  minutes  during 
which  time  the  spacecraft  performs  a  bank  maneuver  to  reduce  the  loads  on  the  vehicle. 
A  plot  of  bank  angle  versus  time  (Figure  49)  shows  this  maneuver  in  greater  detail.  As 
before,  the  bank  angle  begins  lift -up  until  approximately  1.4  minutes  into  the  trajectory. 
At  that  time  the  vehicle  banks  instantaneously  to  lift -down  for  30  seconds  before 
switching  back  to  lift-up  for  another  20  seconds.  Finally,  the  bank  angle  switches  lift- 
down  for  the  remainder  of  the  trajectory. 
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Velocity  History 
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Figure  48:  Velocity  History  (G-limited) 
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Figure  49:  Control  History  (G-limited) 
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The  first  bank  angle  reversal  is  clearly  reducing  the  tangential  loads  on  the 
vehicle,  as  shown  in  the  following  plots  of  body  accelerations.  The  first  lift -down 
segment  corresponds  exactly  to  the  time  interval  for  which  tangential  acceleration  is 
nearly  constant  at  the  constrained  3-g  limit.  This  maneuver  shows  corresponding 
switches  in  the  normal  and  bi-normal  accelerations,  coincident  with  the  bank  angle 
maneuvers.  The  second  bank  angle  correction  is  more  intriguing  as  it  occurs  after  the 
tangential  loads  are  decreasing  in  magnitude.  This  maneuver  seems  to  be  adding  a  slight 
increase  in  flight  path  angle  to  offset  the  period  in  which  the  flight  path  angle  was 
relatively  constant  while  lift -down  (Figure  48). 
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Figure  50:  Acceleration  Histories  (G-limited) 


The  activation  of  the  g-limit  constraint  can  be  further  verified  by  examination  of 
the  covector  associated  with  the  constraint  that  is  provided  by  DIDO.  Figure  51  clearly 
shows  that  the  constraint  becomes  active  over  the  time  period  of  the  first  bank  reversal 
maneuver. 
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Figure  51:  G- limit  Constraint  Covector  History 
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The  heating  rate  is  less  than  the  unconstrained  case,  peaking  at  30.39  J/cnr 
(unconstrained  value  of  35.02  J/cm2)  but  the  total  heat  load  is  higher  at  3814.6  W/cm2 
(unconstrained  value  of  3388.4  W/cm').  This  is  consistent  with  shallower  trajectories 
which  tend  to  have  longer  but  smaller  heat  pulses.  This  higher  heat  load  contributes  to  a 
slightly  higher  performance  index  as  total  heat  load  is  mapped  to  heat  shield  mass  in  our 
model.  The  total  aerocapture  mass  increases  from  28.75  kg  to  30.2  kg  (Table  10). 


Propellant  mass 

11.13  kg 

Front-shield  mass 

19.07  kg 

Total 

30.20  kg 

Table  10:  Cost  Function  Breakdown  (G-limited) 
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F.  MAXIMUM  AEROCAPTURE  CORRIDOR  SUBJECT  TO  HEATING 

RATE  CONSTRAINT 

One  of  the  difficulties  with  implementation  of  the  aerocapture  concept  is  the 
precision  with  which  the  spacecraft  must  guided  to  atmospheric  interface.  Aerocapture 
corridors  are  typically  quite  narrow,  measuring  only  a  couple  degrees.  The  upper  limit 
(that  is  to  say  shallowest  angle)  is  normally  defined  as  the  shallowest  initial  flight  path 
angle  for  which  a  lift -down  bank  angle  profile  will  successfully  meet  the  terminal 
conditions  [Ref  9].  Similarly  the  lower  limit  is  defined  as  the  initial  flight  path  angle  for 
which  a  lift-up  bank  angle  profile  will  meet  the  terminal  conditions.  Some  preliminary 
work  was  done  in  investigating  whether  an  optimal  bank  angle  strategy  can  increase  the 
aerocapture  corridor  width. 

The  maximum  aerocapture  corridor  problem  was  solved  by  separately  solving  two 
related  problems:  maximum  initial  flight  path  angle  and  minimum  initial  flight  path 
angle;  respectively: 


J=~  Yo 

(183) 

J=  Yo 

(184) 

The  difference  between  the  minimum  and  maximum  initial  flight  path  angle  is  then  the 
maximum  corridor  width.  In  addition,  both  solutions  were  subject  to  heating-rate 
constraints.  Note  that  it  is  important  that  the  initial  conditions  for  each  problem  be  the 
same  so  that  the  resultant  trajectories  can  be  fairly  compared.  In  fact,  if  this  is  not  done, 
the  optimal  solution  for  the  two  problems  differ  in  initial  heading  by  180  degrees! 
Instead  the  minimum  initial  flight  path  angle  solution  was  generated  first,  and  its  initial 
condition  was  used  to  constrain  the  initial  maximum  flight  path  angle  case.  Moreover, 
with  no  constraint  on  the  final  orbit  inclination,  the  optimal  solutions  placed  the  vehicle 
in  an  equatorial  orbit  -  not  very  desirable  from  a  scientific  point  of  view.  However,  for 
the  equatorial  case,  the  corridor-defining  optimal  bank  profiles  were  constant  lift -up  or 
constant  lift -down  as  per  the  assumption  in  the  definition.  This  turned  out  to  not  be  the 
case  for  trajectories  with  final  inclinations  other  than  zero. 
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Additional  solutions  were  generated  with  a  constraint  on  the  final  orbit  inclination 
of  the  vehicle.  Recalling  that  upper-case  Greek  letters  represent  the  components  of  the 
velocity  vector  resolved  in  the  inertial  frame,  the  inclination  of  the  vehicle  is  related  to 
the  aircraft  states  by: 

cos  i  —  cos  y  cos  'F  y  (185) 

The  maximum  heat  rate  was  set  at  50  W/cmr  and  the  targeted  inclination  was  arbitrarily 
chosen  to  be  70  degrees  such  that 

cos(70deg)<cosd>/  cos'f//  <cos(70deg)  (186) 

Relevant  parameters  for  the  two  trajectories  that  bound  the  maximum  corridor  are 
presented  in  Table  1 1. 
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Minimum  Initial  FPA 

Maximum  Initial  FPA 

Initial  Flight  Path  Angle 

-10.83  deg 

-8.66  deg 

Initial  V- infinity 

5200  m/s 

5200  m/s 

Initial  Latitude 

17.6  deg 

17.6  deg 

Initial  Heading 

-73.4  deg 

-73.46  deg 

Delta-V 

3743.5  m/s 

3662.4  m/s 

Total  Pass  Time 

5.87  min 

11.82  min 

Minimum  Altitude 

56.32  km 

66.52  km 

Max  Dynamic  Pressure 

3298.44  Pa 

1354.9  Pa 

Max  Total  Acceleration 

5.57  g 

2.29  g 

Max  Heating  Rate 

39.76  W/cm2 

27.67  W/cm2 

Heat  Load 

3138.83  J/c2 

4612.22  J/cnr 

Front-shield  Mass 

15.69  kg 

23.06  kg 

Post-Aerocapture  Propellant 

26.41  kg 

11.91  kg 

Total  Aerocapture  Mass 

42.1  kg 

34.97  kg 

Table  11:  Max  Corridor  Boundary  Comparison 


Figure  52  is  a  plot  of  the  position  state  histories  for  the  two  trajectories  (circles 
represent  the  maximum  flight  path  angle  entry  while  the  plus  symbols  represent  the 
minimum  flight  path  angle  entry).  The  trajectories  begin  at  atmospheric  interface  at  the 
optimal  latitude  of  17.6  degrees.  The  minimum  flight  path  angle  trajectory  is  steep 
(-10.83  deg),  with  only  5.87  minutes  of  total  pass  time  compared  with  11.82  minutes  for 
the  shallow  maximum  flight  path  angle  case  (-8.66  deg).  The  steeper  trajectory  naturally 
penetrates  deeper  into  the  atmosphere  to  an  altitude  of  56.32  km  as  compared  with  66.52 
for  the  shallower  trajectory. 


102 


Position  History 


Figure  52:  Max  Corridor  Position  Histories 


The  velocity  histories  are  given  in  Figure  53.  The  total  delta- V  for  the  two 
trajectories  are  nearly  equal  at  about  3703  m/s  +/ -  40  m/s. 


Velocity  History 


Figure  53:  Max  Corridor  Velocity  Histories 
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The  control  histories  for  the  trajectories  are  given  in  the  following  figure.  CMT 
control  histories  were  plotted  in  lieu  of  the  DIDO  solutions  as  they  better  demonstrated 
the  arctangent  characteristic  of  the  bank  angle  schedule.  Clearly  the  optimal  bank  angle 
profiles  are  not  simply  lift-up  or  lift-down  as  was  the  case  for  equatorial  target  orbits. 
The  steep  flight  path  begins  lift -up  and  modulates  to  approximately  60  degrees  during  the 
ascent.  In  a  like  manner,  the  shallower  entry  begins  lift-down  and  modulates  to 
approximately  94  deg  during  the  ascent.  In  neither  case  does  the  lift  vector  cross  the 
horizontal  plane;  instead  the  lift  vector  seems  to  be  used  to  control  depth  of  penetration 
prior  to  the  minimum  altitude  before  switching  to  control  heading  (and  hence  target 
inclination)  for  the  ascent. 
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Figure  54:  Max  Corridor  Control  Histories 


As  expected,  the  peak  heating  rate  for  the  steep  entry  is  significantly  higher  than 
the  shallow  entry  (39.76  W/cm"  as  compared  to  27.67  W/cnr)  although  the  shallower 
trajectory  has  the  higher  heat  load  (4,612  J/cm2  compared  with  3,138  J/cm2)  (Figure  55). 
Note  that  the  peak  heating  rate  was  below  the  constraint  value  of  50  W/cm" 
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Stagnation  Point  Heating 


Figure  55:  Max  Corridor  Heating  Rates 
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Figure  56:  Max  Corridor  Dynamic  Pressures 
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The  peak  dynamic  pressure  (Figure  56)  of  the  steep  trajectory  is  more  than  double 
that  of  the  shallow  trajectory  (3,298  Pa  compared  with  1,355  Pa)  which  leads  to  a  similar 
disparity  in  peak  total  accelerations  (Figure  57).  The  steep  entry  encounters  a  crushing 
5.57  g  peak  acceleration  whereas  the  shallow  entry  peaks  at  a  more  manageable  2.29  g. 
Figure  58  resolves  the  aerodynamic  accelerations  into  flight -path  related  components. 


Acceleration  Mag  vs  Time 


Figure  57:  Max  Corridor  Total  Accelerations 
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Body  Accelerations 
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Figure  58:  Max  Corridor  Acceleration  Components 


The  significant  differences  in  characteristics  of  the  trajectories  which  define  the 
boundaries  of  the  maximum  aerocapture  corridor  illustrate  the  difficulties  imposed  upon 
the  design  team.  To  utilize  the  entire  available  corridor,  the  vehicle’s  structure  must  be 
sized  to  the  more  dynamic,  steep  boundary  trajectory  while  the  TPS  mass  must  be  sized 
for  the  larger  heat  loads  of  the  shallow  entry.  Additionally,  the  TPS  material  selection 
will  be  dependent  on  the  maximum  heating  rate  of  the  steep  trajectory. 
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IX.  FORMULATING  THE  COMBINED  LOW  THRUST  AND 
AEROCAPTURE  PROBLEM 


Previous  chapters  detail  the  problem  formulation  for  both  low  thrust  and 
aerocapture  trajectories  taken  separately.  The  attention  is  now  turned  to  the  problem  of 
solving  for  the  low  thrust  trajectory  that  terminates  with  an  aerocapture  maneuver.  One 
approach  for  generating  feasible  integrated  trajectories  is  to  optimize  the  low  thrust 
trajectory  using  its  final  state  to  derive  the  initial  state  for  an  optimal  aerocapture  pass. 
However,  this  formulation  is  not  optimizing  the  problem  from  end  to  end.  To  find  the 
true  integrated  optimal  solution,  we  must  simultaneously  solve  both  trajectories. 


A.  JUNCTION  CONDITIONS 


Recall  that  a  direct,  discreet  solution  method  reduces  an  optimal  control  problem 
to  a  series  of  dynamic  constraints  sampled  at  various  nodes  (times).  Thus,  the  state 
history  for  a  low  thrust  problem  can  be  expressed  as  a  matrix  such  as 
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Where  the  columns  represent  the  values  of  the  states  at  each  node  time  and  were  j 
is  the  index  corresponding  to  the  final  node  of  the  low  thrust  solution.  Similarly  for 
aerocapture 
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(188) 


where  k  represents  the  final  node  index  of  the  trajectory. 


If  these  two  state  history  matrices  were  of  the  same  dimension,  they  could  be 
augmented  forming  a  combined  state  history  for  the  entire  problem.  However,  note  that 
the  low  thrust  state  history  matrix  has  one  less  row  (state)  than  the  aerocapture  matrix. 
This  can  be  resolved  by  simply  adding  a  “dummy  variable”  to  the  last  row  of  the  matrix. 
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Our  combined  state  history  matrix  now  takes  the  desired  form  of  a  6  by  §+k ) 
matrix  with  the  final  column  of  the  low  thrust  trajectory  occurring  at  index  j  and  the  first 
column  of  the  aerocapture  trajectory  beginning  with  column  /+/. 
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However,  we  now  must  somehow  match  the  physical  meaning  of  the  end  of  the 
low  thrust  trajectory  with  that  of  the  beginning  of  the  aerocapture  trajectory.  As 
formulated,  the  states  take  on  very  different  meanings  in  the  two  portions  of  the  problem. 
For  example,  the  5h  row  of  the  low  thrust  portion  represents  the  vehicle’s  mass  while  in 
the  aerocapture  portion  that  same  row  represents  heading  angle.  Thus  we  need  some 
junction  conditions  to  relate  the  variables  on  the  two  sides  of  the  problem. 
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First,  recall  from  previous  chapters  that  the  initial  final  radius  of  the  low  thrust 
portion  and  the  initial  radius  of  the  aerocapture  problem  are  fixed.  That  is 


r(tj)  =  r 


planet 


T 

atmos  limit 


(191) 

(192) 


Final  velocities  for  the  low  thrust  portion  are  left  unconstrained  to  allow  for  non¬ 
rendezvous  (excess  arrival  V-infinity)  trajectories  to  be  generated.  The  vehicle  mass  at 
the  end  of  the  low  thrust  portion  remains  subject  to  the  physical  limitation  that 
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Although  mass  is  not  a  state  for  the  aerocapture  segment,  its  value  is  needed  for 
the  dynamics  calculations.  Assuming  that  there  is  no  post  low  thrust  staging  this  yields 


m  ,  =  m , 

aerocapture  \  j 


(194) 


Since  vehicle  mass  at  arrival  is  on  the  same  order  as  the  initial  vehicle  mass  at 
launch,  there  conveniently  is  no  need  rescale  the  problem  as  the  normalizing  units  of 
mass  may  be  chosen  to  be  the  same  for  both  trajectory  segments. 

Recall  that  for  convenience  the  initial  longitude  of  the  aerocapture  trajectory  was 
set  to  zero  so 
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Similarly  the  initial  velocity  states  of  the  aerocapture  problem  can  be  related  in 
some  manner  to  the  terminal  conditions  of  the  low  thrust  problem. 
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(198) 
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To  see  exactly  how  these  are  related,  the  V-infinity  of  arrival  for  the  aerocapture 
problem  must  be  found.  The  heliocentric  inertial  velocity  of  the  vehicle  is  give  by  the 
vector  equation 
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which  can  be  rearranged  to  find  the  inertial  velocity  of  the  vehicle  relative  to  the  planer 
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With  respect  to  the  variables  used  in  the  low  thrust  portion,  the  velocity  vector  of 
the  vehicle  with  respect  to  the  sun  is 


and  the  circular  velocity  of  the  planet  is  given  by 
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Thus  the  velocity  of  the  vehicle  relative  to  the  planet  can  be  expressed  as 
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The  magnitude  of  this  vector  is 
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The  square  root  term  leads  to  possibility  of  a  singularity  so  it  is  more  convenient  to  use 
the  square  of  the  V-infinity  at  arrival  given  as 
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At  this  point,  the  inertial  arrival  velocity  magnitude  in  the  planet  frame  at  the 
beginning  of  the  aerocapture  trajectory  is  known.  Given  the  small  size  of  the  target 
planet  relative  to  the  scale  of  the  heliocentric  trajectory,  small  course  corrections  far  from 
the  target  planet  allow  for  any  point  on  the  B-plane  to  be  targeted  (the  B-plane  is  a 
reference  plane  used  for  interplanetary  targeting).  This  means  that  the  initial  latitude, 
longitude,  velocity,  heading  and  flight  path  angle  sates  for  the  aerocapture  problem  is  are 
essentially  free,  provided  they  all  take  values  that  are  consistent  with  the  arrival  V- 
infinity. 

Using  vis -viva,  the  magnitude  of  the  velocity  vector  at  the  atmospheric  interface 
(where  we  begin  the  aerocapture  trajectory)  is  related  to  the  arrival  V- infinity. 

- f^  +  JL  =  o  (206) 
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Again  referring  to  Appendix  A,  the  velocity  components  in  a  rotating  REN  frame 
may  be  related  to  the  arrival  V-infinity,  thus  we  can  now  relate  the  V-infinity  as 
calculated  from  the  final  states  of  the  low  thrust  trajectory  to  the  initial  states  of  the 
aerocapture  trajectory  as 

M|“  =v(^+1)2+r(r,+i)^2cos2^(r.+1)---  (2Q7) 

+  2r  (tj+l )  v  (tj+1 )  f2cos(^  (tj+1 )  cos v  (tJ+l )  easy  (t j+l ) 


Substituting  Eqn.  (205)  into  the  above  equation  yields  the  following  important  junction 
condition 
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However,  recall  that  not  only  were  the  states  different  between  the  two 
formulations,  but  they  were  scaled  markedly  different  as  well.  This  can  be  resolved 
through  the  addition  of  a  constant  conversion  factor 
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+  2 r  (tj+l )  v  (tj+x )  Q.  cos(^  (tj+l )  cosy  (tJ+ , )  cos y  (f;+1 ) 


Eqns.  ( 191)-(  193),  (195),  and  (209)  complete  the  junction  conditions  for  the 
combined  problem. 

B.  COST  FUNCTIONS 

1.  Minimum  Total  Propellant 

A  relatively  simple  to  implement  performance  index  for  the  combined  problem  is 
to  minimize  the  total  propellant  required  for  the  mission.  This  can  be  accomplished 
despite  the  propellants  being  of  different  types  with  different  Isp.  This  combined 
propellant  cost  can  be  expressed  in  Mayer  form  as 

J  =  —  (m  .  ,  +m  .  )  (210) 

\  proplowthrust  prop  circ  /  v  7 

where  m proplowthrust  is  the  propellant  mass  consumed  during  the  low  thrust  portion  of  the 
trajectory  as  given  by 

^ proplowthrust  ^  ^ j  )  "•(<„)  (211) 

and  mpropcirc  is  the  mass  of  propellant  needed  to  circularize  the  post-aerocapture  transfer 
orbit. 


2.  Maximum  Scientific  Mass 

The  total  maximum  scientific  mass  at  arrival  can  found  by  taking  the  total 
propellant  cost  above  and  adding  to  it  the  estimated  mass  of  the  heat  shield.  Using  the 
same  mapping  between  heat  load  and  front -shield  mass  used  in  chapters  VII  and  VIII  this 
can  be  expressed  as  a  Bolza  cost  as  follows 
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X.  SOLUTIONS  TO  THE  COM  BINED  LOW  THRUST  AND 
AEROCAPTURE  PROBLEM 


A.  MINIMUM  TOTAL  FUEL  SOLUTION 

The  minimum  total  fuel  solution  was  obtained  using  40  nodes  of  resolution  for 
each  segment  of  the  trajectory.  Obtaining  feasible  trajectories  proved  rather  difficult. 
Increasing  the  thrust  capacity  of  the  low  thrust  vehicle  helped  to  obtain  solutions  at  the 
expense  of  realism.  In  addition,  the  rotation  of  the  target  planet  was  set  to  zero  to  simplify 
the  solution  and  reduce  computation  time.  This  was  accomplished  by  simply  assuming 
that  the  vehicle  employed  10  NSTAR  engines.  As  show  in  Figure  59  (where  again  the 
circles  are  the  DIDO  solution  and  the  solid  line  is  the  propagated  trajectory)  the  trajectory 
begins  with  zero  Q  and  a  700  kg  vehicle  and  arrives  at  1.52  AU  291  days  after  launch. 
In  the  process,  61.58  kg  of  propellant  are  consumed  for  an  arrival  mass  at  Mars  of 
638.42  kg.  The  close  agreement  between  the  DIDO  solution  and  the  propagated  solution 
verify  the  feasibility  of  the  solution. 
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Low  Thrust  States:  Earth  to  Mars 


Figure  59:  Low  Thrust  State  Histories  (Min  Total  Fuel) 


The  heliocentric  trajectory  is  shown  in  Figure  60.  An  initial  burn  of 
approximately  48  days  increases  the  semi -major  axis  of  the  transfer  orbit  until  the 
aphelion  intersects  Mars’s  orbit.  The  small  deviations  between  the  DIDO  and  propagated 
solutions  are  more  visible  in  this  presentation. 
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Low  Thrust:  Earth  to  Mars 


Figure  60:  Heliocentric  Low  Thrust  Trajectory  (Min  Total  Fuel) 


The  vehicle  arrives  at  Mars  with  an  inertial  arrival  V-infinity  of  2626.2  m/s  and  a 
relatively  shallow  flight  path  angle  of  -7.62  degrees.  Note  that  this  is  a  relatively  shallow 
flight  path  angle  but  the  correspondingly  small  arrival  Vinfinity  ensures  that  the  vehicle 
will  not  “skip”  off  the  atmosphere.  Figure  61  and  Figure  62  depict  the  vehicles  states  for 
the  aerocapture  portion  of  the  trajectory.  A  minimum  altitude  of  70  km  is  reached  3 
minutes  into  the  13.65  minute  trajectory.  The  pass  yields  a  total  delta-V  of  2078  m/s. 
Again  the  DIDO  solution  and  the  propagated  solution  are  in  excellent  agreement.  The 
transfer  orbit  is  highly  circular  with  an  eccentricity  of  only  0.03.  This  contributes  to  the 
low  propellant  mass  required  to  circularize  of  only  12.1  kg. 
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Figure  61:  Aerocapture  Position  States  (Min  Total  Fuel) 

Velocity  History 
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Figure  62:  Aerocapture  Velocity  States  (Min  Total  Fuel) 
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The  bank  angle  history  is  given  in  Figure  63.  Unfortunately,  current  versions  of 
DIDO  do  not  return  the  covectors  for  problems  with  internal  knots  such  as  this  one 
making  verification  difficult.  For  this  reason,  the  CMT  controls  can  not  be  shown  for 
comparison.  Similarly  the  Hamiltonian  is  not  available  for  verification  of  the  first 
integral. 
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Figure  63:  Aerocapture  Control  History  (Min  Total  Fuel) 


Figure  64  gives  the  heating  rate  over  the  trajectory.  The  peak  heating  rate  of  12. 2 

9  9 

W/cnr  occurs  2.14  minutes  into  the  trajectory  and  the  total  heat  load  is  3054.3  J/cnr. 
Using  50  kg  of  TPS  mass  per  10,000  J/cm2  of  heat  load,  this  corresponds  to 
approximately  15.3  kg  of  front-shield  mass. 
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Stagnation  Point  Heating 
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Figure  64:  Aerocapture  Heating  Rate  (Min  Total  Fuel) 


Unfortunately,  DIDO  does  not  currently  return  covectors  for  problems  with 
interior  knots  such  as  this  mixed -dynamics  problem.  Without  the  co  vectors  or 
Hamiltonian,  verification  of  the  optimality  of  the  solution  is  more  difficult.  In  principle, 
the  interior  event  conditions  of  the  mixed-dynamic  problem  could  be  used  to  formulate 
two  optimization  problems,  essentially  breaking  the  problem  back  into  its  parts.  Each 
optimal  control  problem  could  then  be  solved,  and  the  covectors  and  Hamiltonians 
exploited  to  determine  optimality.  The  optimality  of  the  combined  problem  could  then  be 
declared  if  each  individual  problem  was  optimal  on  its  own  and  each  solution  state  and 
control  history  matched  that  of  the  combined  solution.  This  verification  was  not 
performed  in  this  work  due  to  lack  of  time. 
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B.  MAXIMUM  SCIENTIFIC  MASS 


The  maximum  scientific  mass  cost  function  added  further  difficulties  in  obtaining 
a  solution.  Simply  changing  the  cost  function  yielded  infeasible  solutions.  Unlike  the 
minimum  total  fuel  case,  increasing  the  vehicle  thrust  did  not  resolve  this  issue.  To 
obtain  a  solution,  a  hypothetical  target  planet  with  the  properties  of  Mars  was  placed  at 
5.2  AU.  Like  the  minimum  total  fuel  case,  the  problem  was  solved  for  a  non-rotating 
target  atmosphere.  Even  with  these  changes,  solutions  for  high  numbers  of  nodes  became 
numerically  unstable  and  yielded  infeasible  solutions.  The  solution  presented  below  was 
obtained  using  120  nodes  for  the  low  thrust  segment  and  50  nodes  for  the  aerocapture 
segment.  The  state  histories  for  this  solution  are  given  in  Figure  65.  The  vehicle  begins 
with  a  zero  C3  and  initial  mass  of  660  kg  and  arrives  at  Mars  with  487  kg  consuming  173 
kg  of  propellant.  The  total  trip  time  is  1187.4  days.  The  DIDO  low  thrust  state  histories 
are  in  excellent  agreement  with  the  propagated  states. 


Low  Thrust  States:  Earth  to  5  AU 


Figure  65:  Low  Thrust  State  Histories  (Maximum  Scientific  Mass) 
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Figure  66  shows  the  heliocentric  trajectory  for  the  low  thrust  portion.  Again,  the 
vehicle  conducts  a  long  continuous  burn  for  175  days  before  shutting  off  and  coasting 
with  just  enough  energy  to  reach  5.2  AU. 


Low  Thrust:  Earth  to  5  AU 


AU 

Figure  66:  Heliocentric  Low  Thrust  Transfer  Orbit  (Max  Scientific  Mass) 


Figure  67  shows  the  control  histories  for  the  low  thrust  segment  of  the  mission. 
The  lower  plot  shows  the  normalized  thrust  history  with  distinct  thrust  switch  at 
approximately  175  days.  The  thrust  angle  is  approximately  zero  during  the  thrusting 
portion  of  the  trajectory.  The  remainder  of  the  thrust  angle  history  may  be  disregarded  as 
it  is  of  course  meaningless  when  thrust  magnitude  is  zero. 
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Figure  67:  Low  Thrust  Control  Histories  (Max  Scientific  Mass) 


The  aerocapture  state  histories  are  show  in  Figure  68  and  Figure  69.  The 
trajectory  begins  with  an  inertial  Y-infinity  at  arrival  of  5211  m/s  and  a  flight  path  angle 
of  arrival  of  -10.1  degrees.  The  fact  that  the  initial  Yinfinity  is  not  zero  demonstrates 
that  the  global  solution  is  minimizing  low  thrust  propellant  at  the  expense  of  more 
efficient  thermal  energy  dissipation  during  the  aerocapture  segment.  The  initial  heading 
of  -95.6  degrees  is  due  to  the  non-rotating  atmosphere  which  causes  the  cost  function  to 
be  invariant  with  latitude  and  heading  angle.  The  inertial  delta- V  for  the  pass  is  3668  m/s 
and  terminates  in  an  orbit  with  an  eccentricity  of  0.04.  Capture  occurs  at  1.87  minutes. 
Unlike  the  low  thrust  trajectories,  the  DIDO  solutions  do  not  correspond  as  well  with  the 
propagated  states,  particularly  with  the  radius  state.  It  is  likely  that  increasing  the  number 
of  nodes  for  this  segment  would  lead  to  better  convergence  between  the  DIDO  and 
propagated  solutions. 
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Figure  68:  Aerocapture  Position  States  (Max  Scientific  Mass) 
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Figure  69:  Aerocapture  Velocity  States  (Max  Scientific  Mass) 
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The  bank  angle  control  history  is  shown  in  Figure  70  below.  As  in  the  pure 
aerocapture  optimizations  investigated  earlier,  the  bank  angle  begins  approximately  lift 
up  before  switching  at  approximately  1.9  minutes  to  lift  down.  Again,  the  chatter  in  the 
early  history  of  the  control  is  due  to  the  lack  of  aerodynamic  control  authority  high  in  the 
atmosphere  and  would  be  expected  to  be  smoothed  with  higher-node  solutions. 


Control  History 


Figure  70:  Bank  Angle  History  (Max  Scientific  Mass) 


Figure  71  shows  the  heating  rate  history  during  the  atmospheric  pass.  The 
maximum  stagnation  point  heating  of  34.2  W/cm2  occurs  at  1.24  minutes  and  the  total 
heat  load  is  3586.59  J/cm2.  Using  the  TPS  mass  model  discussed  previously,  this 
corresponds  to  a  front-shield  mass  of  17.93  kg.  Combined  with  the  9.83  kg  of  propellant 
required  to  circularize  the  orbit,  the  total  aerocapture  mass  is  27.7  kg.  This  compares 
with  331  kg  that  would  be  required  for  a  pure  propulsive  capture  at  the  target  planet. 
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Figure  71:  Aerocapture  Heating  Rate  (Max  Scientific  Mass) 
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XL  CONCLUSIONS 


The  suitability  of  the  direct  pseudospectral  method  for  solving  highly  non-linear 
astrodynamic  problems  has  been  explored.  For  low  thrust  problems,  the  method  has  been 
shown  to  produce  excellent  results,  particularly  for  minimum  time  problems.  However 
unknown  factors  cause  low  thrust  minimum  fuel  solutions  to  be  more  difficult  to 
consistently  obtain.  The  method  was  particularly  successful  in  solving  optimal 
aerocapture  trajectories  over  a  range  of  conditions. 

The  suitability  of  the  direct  method  for  simultaneously  solving  a  combined  low 
thrust  trajectory  with  terminal  aerocapture  was  also  explored.  Although  the  fidelity  of  the 
models  was  reduced  to  obtain  feasible  solutions,  the  proof-of-concept  is  considered  a 
success  as  it  successfully  found  feasible  solutions  for  the  combined  trajectories.  This 
concept  of  simultaneously  optimizing  trajectory  segments  with  vastly  different  dynamics 
has  the  potential  to  identify  previously  unexplored  trajectory  combinations  and  further 
research  in  this  area  is  strongly  suggested. 
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XH.  FUTURE  WORK 


Due  to  the  proof-of-concept  nature  of  this  work,  there  are  numerous  areas  for 
which  future  work  is  encouraged.  For  low  thrust  problems,  launch  vehicle  optimization 
can  be  added  allowing  for  design  trades  between  initial  mass  and  initial  C3.  More 
realistic  trajectories  may  be  obtained  by  exploring  3-DOF  as  well  as  non-circular  target 
orbits.  Additionally,  adding  the  gravitational  effects  of  additional  bodies  may  allow  for 
the  exploitation  of  gravity  assists  and  further  mass  savings.  Finally,  the  difficulties  with 
obtaining  certain  low  thrust  minimum  fuel  trajectories  should  be  further  investigated. 

For  aerocapture,  a  more  accurate  TPS  mass  model  should  be  developed  such  that 
both  the  heat  load  and  the  peak  heating  rate  are  taken  into  account.  Furthermore,  the 
benefits  of  angle  of  attack  modulation  during  aerocapture  as  well  as  thrusting  arcs  may 
yield  new  families  of  trajectories  and  should  be  explored.  Additionally,  other  cost 
functions  such  as  minimum  altitude  may  be  useful  for  such  missions  as  sample  return. 

Finally,  this  initial  work  solving  mixed  dynamic  optimization  problems  may  be 
expanded  in  many  areas.  To  obtain  any  solutions  at  all,  the  fidelity  was  considerably 
reduced.  Hopefully  many  of  the  problems  with  the  combined  trajectories  will  be  rectified 
when  the  difficulties  with  the  minimum  fuel  trajectories  discussed  above  are  resolved. 
Lastly,  with  future  versions  of  DIDO,  users  should  be  able  to  recover  the  covectors  for 
problems  with  interior  knots,  allowing  for  the  verification  of  the  DIDO  trajectories. 
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APPENDIX  A:  USEFUL  TRANSFORMATIONS 


A.  COORDINATE  TRANSFORMATIONS 

The  following  transformations  are  useful  for  moving  between  various  aerocapture 

frames. 

1.  Spherical  to  Inertial 

First  rotate  about  K  by  0  ,  then  Y '  by  -<|) . 


er 

COS(|) 

0 

—  sin  ([) 

COS0 

sin0 

O' 

= 

0 

1 

0 

-sin0 

COS0 

0 

J 

e. 

L  <t>  J 

sintf) 

0 

cost]) 

0 

0 

1 

K 

which  simplifies  to 


COS0  cose]) 

sin0  cost]) 

sincf) 

-sin0 

COS0 

0 

J 

-cos0  sincf) 

-sin0  sin  ([) 

cost]) 

K 

2.  REN  Frame  to  Frenet  Frame 

First  rotate  about  r  by\|/  ,  then  rotate  about  cf) '  by  — y  . 


en 

cosy 

-siny 

O' 

'1 

0 

0 

= 

siny 

cosy 

0 

0 

cos))/ 

sin))/ 

_K_ 

0 

0 

1 

0 

-sin\)/ 

COS))/ 

e. 

L  $  J 

which  simp  lifies  to 

cosy 
siny 
0 


-siny  cos))/ 

-siny  sin))/ 

cosy  cos))/ 

cosy  sin))/ 

-sin))/ 

COS))/ 

e. 

L  <t>  J 

which  can  be  inverted  to 


cosy 

siny 

0 

-siny  cos))/ 

cosy  cos))/ 

-sin))/ 

-siny  sin))/ 

cosy  sin))/ 

cos))/ 

_K_ 

(213) 


(214) 


(215) 


(216) 


(217) 
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which  can  be  inverted  to 


COS0  cost]) 
sin0  cost]) 
sin  tf) 


-sin0 

-cos0  sintf) 

COS0 

-sin0  sintf) 

% 

0 

cost]) 

e. 

L  <l>  J 

(218) 


3.  PCF  to  PCI  Velocity 

Recall  that  Ypa  is  the  inertial  velocity  vector  and  Ypcf  is  the  velocity  vector  in 
the  rotating  frame.  The  transport  theorem  states 


V  —V  +PCI  9  PCF  xr 
v  pci  v  pcf  +  •  xr 


But  the  angular  rate  of  the  rotating  frame  is  given  by 


PCI  ?  PLF  =  Q.K  =  sin  tf)  er  +  cost])  \ 


The  positions  vector  in  spherical  coordinates  is  simply 

{ r]  =  r  <?„ 

The  can  be  transformed  from  REN  to  spherical  by 

{VpcU  = 


cosy 

siny 

0 

"0" 

-  sin  y  cost)/ 

cosy  cost)/ 

-sint)/ 

V 

-siny  sint)/ 

cosy  sint)/ 

cost)/ 

_0_ 

which  gives 


|  V pep }  =  vsiny  er  +  vcosy  cost)/  ee  +v  cosy  sint)/  e ^ 


so  substituting  into  Eqn.  (219)  results  in 

{VpC/Irftj, 


vsiny 

e 

r 

4 

e, 

4> 

vcosy  cost)/ 

+ 

Q  sintf) 

0 

£2cOSt|) 

vcosy  sint)/ 

r 

0 

0 

(219) 


(220) 


(221) 


(222) 


(223) 


(224) 


which  simplifies  to 

=  vsiny  er  +  (vcosy  cost)/  +/T2  cost]))  e6  +v  cosy  sint)/  e ^  (225) 


134 


This  gives  the  spherical  inertial  velocity  components  Vr ,  Ve  ,  and  Vti)  in  terms  of 
the  velocity  vector  in  the  rotating  REN  frame. 

R}ra  =vsiny 

{V9 } />a  =v cosy  cost)/  +r£2c os(f)  (226) 

K}ra  =v cosy  sin it)/ 


We  can  now  obtain  the  velocity  components  in  the  inertial  REN  frame  by  squaring  both 
sides  of  Eqn.  (225)  to  get 

||ypc/|r  =  v2  +  r'Qr  cos2  cj)  +  2rvf2 cos (J)  cost)/  cosy  (227) 

so  the  inertial  speed  in  the  REN  frame  is 

VPCI  =  yjv2  +  r2Cl2  c os2(|)  +2  rv'f2cos(|)  cosy  cosy  (228) 


Recalling  Figure  3 


VE  =  atan 


(229) 


and 


f 


r  =  atan 


V 


(230) 


Substituting  Eqns.  (226)  into  (229)  and  (230)  we  get 


( 


T*  =  atan 


vcosy  sin\)/ 


v cosy  cost)/  +r£2cos(J> 


(231) 


and 


r  =  atan 


vsiny 


Jcos2y  (v2  +fl2r2)+  2v/i2cosy  cost)/  cos  (|) 


(232) 


135 


Thus  Eqns  (228),  (231),  and  (232)  give  us  the  inertial  velocity  components 
resolved  in  the  REN  frame  as  functions  of  the  position  and  velocity  components  as 
measured  in  the  rotating  REN  frame. 
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APPENDIX  B:  MINIMUM  AEROCAPTURE  MASS  AT  MARS 
FROM  ZERO  ARRIVAL  V-INFINITY  DATA 


Data  Summary: 

Inertial  Velocity  Components: 

Arrival  V  infinity  (m/s):  0.00 

Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Delta-V  (m/ s)  : 

Rotating  Velocity  Components 
Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 


Atmospheric  Exit: 

Speed  (m/ s) : 

3274.53 

Heading  (deg) : 

0.09 

Flight  Path  Angle  (deg) : 

2  .  02 

Delta-V  (m/ s)  : 

1427 . 83 

4702.36 

-0.09 

-7.83 


3523.53 

0.08 

1.87 

1425.77 


4949.30 

-0.09 

-7.44 


Trajectory  Analysis: 
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Total  Pass  Time: 


734.40  sec  (12.24  min) 


Minimum  Altitude  (km) : 

Time  to  Min  Alt : 

Max  Dynamic  Pressure  (Pa) : 

Time  of  Max  Dynamic  Pressure: 
Max  Acceleration:  (g) 

Time  of  Max  Acceleration: 

Max  Heating  Rate:  (W/cmA2) 

Time  of  Max  Heating  Rate: 

Heat  Load:  (J/cmA2) 

Final  Orbit  Parameters: 
Semi-major  Axis  (km) : 

Periapsis  (km) : 

Apoapsis  (km) : 

Apoapsis  Altitude  (km) : 
Eccentricity: 

Capture  Time: 

Total  Aerocapture  Mass  (kg) : 
Propellant  Mass  to  Circularize 
Estimated  Front-shield  mass 


70 . 14 

158.94  sec  (2.65  min) 

505.71 

158.94  sec  (2.65  min) 

0.85 

158.94  sec  (2.65  min) 

7.44 

140.18  sec  (2.34  min) 

1700.38 

3563.65 

3437.39 
3689.92 
300.00 
0.04 

0.00  sec  (0.00  min) 

19.14 

After  Aerocapture  (kg):  10.63 

) :  8.50 


Pure  Propulsive  insertion: 

Isp  =  330.000000 

Propellant  Mass  (kg):  201.442991 
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APPENDIX  C:  MINIMUM  AEROCAPTURE  MASS  AT  MARS 
FROM  EXCESS  ARRIVAL  V-INFINITY  DATA 


Data  Summary: 

Inertial  Velocity  Components: 

Arrival  V  infinity  (m/s) : 

Atmospheric  Entry: 

Speed  (m/s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 
Delta-V  (m/ s) : 

Rotating  Velocity  Components: 
Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 
Delta-V  (m/ s)  : 


5200.00 

7178 . 83 

-0.15 

-10.47 

3517.05 
0.27 
2 . 18 
3661.78 


6933.98 

-0.15 

-10.85 

3268.11 

0.29 

2.35 

3665.88 


Trajectory  Analysis: 


139 


Total  Pass  Time: 

Minimum  Altitude  (km) : 

Time  to  Min  Alt: 

Max  Dynamic  Pressure  (Pa) : 

Time  of  Max  Dynamic  Pressure: 
Max  Acceleration:  (g) 

Time  of  Max  Acceleration: 

Max  Heating  Rate:  (W/cmA2) 

Time  of  Max  Heating  Rate: 

Heat  Load:  (J/cmA2) 

Final  Orbit  Parameters: 
Semi-major  Axis  (km) : 

Periapsis  (km) : 

Apoapsis  (km) : 

Apoapsis  Altitude  (km)  : 
Eccentricity : 

Capture  Time: 

Total  Aerocapture  Mass  (kg) : 
Propellant  Mass  to  Circularize 
Estimated  Front-shield  mass  (k< 


601.02  sec  (10.02  min) 

58.51 

94.63  sec  (1.58  min) 

2703.70 

79.78  sec  (1.33  min) 

4.57 

79.78  sec  (1.33  min) 

35.02 

72.76  sec  (1.21  min) 

3388.37 

3550.25 

3410.57 
3689.92 
300.00 
0.04 

110.50  sec  (1.84  min) 

28.75 

After  Aerocapture  (kg) :  11.81 

):  16.94 


Pure  Propulsive  insertion: 

Isp  =  330.000000 

Propellant  Mass  (kg):  386.401888 
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APPENDIX  D:  MINIMUM  AEROCAPTURE  MASS  AT  MARS 
FROM  EXCESS  ARRIVAL  V-INFINITY  DATA 


Data  Summary: 

Inertial  Velocity  Components: 

Arrival  V  infinity  (m/s) : 

Atmospheric  Entry: 

Speed  (m/s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 

Atmospheric  Exit: 

Speed  (m/s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 
Delta-V  (m/ s)  : 

Rotating  Velocity  Components: 
Atmospheric  Entry: 

Speed  (m/s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 
Delta-V  (m/ s)  : 


9419.  99 

24987.43 

-0.08 

-9.35 

16143.37 

0.54 

1.28 

8844.06 


22615.43 

-0.08 

-10.34 

13736.80 

0.64 

1.51 

8878.62 


Trajectory  Analysis: 
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Total  Pass  Time: 

1358.47  sec 

(22. 

64  min) 

Minimum  Altitude  (km) : 

256.46 

Time  to  Min  Alt: 

249.75  sec 

(4.16 

min) 

Max  Dynamic  Pressure  (Pa) : 

3451.29 

Time  of  Max  Dynamic  Pressure: 

231.55  sec 

(3.86 

min) 

Max  Acceleration:  (g) 

5.83 

Time  of  Max  Acceleration: 

231.55  sec 

(3.86 

min) 

Max  Heating  Rate:  (W/cmA2) 

264.59 

Time  of  Max  Heating  Rate: 

213.90  sec 

(3.56 

min) 

Heat  Load:  (J/cmA2) 

55903.98 

Final  Orbit  Parameters: 

Semi-major  Axis  (km) : 

24757.55 

Periapsis  (km) : 

23891 . 11 

Apoapsis  (km) : 

25624.00 

Apoapsis  Altitude  (km) : 

1000.00 

Eccentricity : 

0.03 

Capture  Time: 

213.90  sec 

(3.56 

min) 

327.98 
48 . 40 
279.52 

Pure  Propulsive  insertion: 

Isp  =  330.000000 

Propellant  Mass  (kg):  528.735021 


Total  Aerocapture  Mass  (kg) : 

Propellant  Mass  to  Circularize  After  Aerocapture  (kg) : 
Estimated  Front-shield  mass  (kg) : 
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APPENDIX  D:  MINIMUM  AEROCAPTURE  MASS  AT  MARS 

WITH  G- LIMITS 

Data  Summary: 

Inertial  Velocity  Components: 

Arrival  V  infinity  (m/s):  5199.99 

Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg)  : 

Flight  Path  Angle  (deg) 

Delta-V  (m/ s)  : 

Rotating  Velocity  Components 
Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 


Atmospheric  Exit: 

Speed  (m/ s) : 

3271 . 83 

Heading  (deg) : 

-0.32 

Flight  Path  Angle  (deg) : 

2.16 

Delta-V  (m/ s)  : 

3661.56 

Trajectory  Analysis: 

6933.39 

-0.09 

-10.07 


3520.80 

-0.29 

2.01 

3658.02 


7178 . 82 

-0.09 

-9.73 


Total  Pass  Time: 


684.58  sec  (11.41  min) 
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Minimum  Altitude  (km) : 

Time  to  Min  Alt: 

Max  Dynamic  Pressure  (Pa) : 

Time  of  Max  Dynamic  Pressure: 
Max  Acceleration:  (g) 

Time  of  Max  Acceleration: 

Max  Heating  Rate:  (W/cmA2) 

Time  of  Max  Heating  Rate: 

Heat  Load:  (J/cmA2) 

Final  Orbit  Parameters: 
Semi-major  Axis  (km) : 

Periapsis  (km) : 

Apoapsis  (km) : 

Apoapsis  Altitude  (km)  : 
Eccentricity : 

Capture  Time: 

Total  Aerocapture  Mass  (kg) : 
Propellant  Mass  to  Circularize 
Estimated  Front-shield  mass  (k< 


61.98 

122.65  sec  (2.04  min) 

1803.97 

96.30  sec  (1.61  min) 

3.05 

96.30  sec  (1.61  min) 

30.39 

78.32  sec  (1.31  min) 

3814.61 

3557.99 

3426.09 

3689.90 

299.98 
0.04 

136.75  sec  (2.28  min) 

30.20 

After  Aerocapture  (kg):  11.13 

r)  :  19.07 


Pure  Propulsive  insertion: 

Isp  =  330.000000 

Propellant  Mass  (kg):  386.401510 
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APPENDIX  E:  MAXIMUM  INITIAL  FLIGHT  PATH  ANGLE  FOR 

AEROCAPTURE  AT  MARS 

Data  Summary: 

Inertial  Velocity  Components: 

Arrival  V  infinity  (m/s) :  5200.00 

Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Atmospheric  Exit: 

Speed  (m/s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Delta-V  (m/ s)  : 

Rotating  Velocity  Components 
Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Delta-V  (m/ s)  : 

Trajectory  Analysis: 

Total  Pass  Time:  708.99  sec  (11.82  min) 


3436.81 

-69.58 

2.26 

3678.72 


7115.54 

-75.32 

-8.73 


3516.46 

-66.33 

2.21 

3662.37 


7178 . 83 

-73.46 

-8.65 
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Minimum  Altitude  (km) : 

Time  to  Min  Alt: 

Max  Dynamic  Pressure  (Pa) : 
Time  of  Max  Dynamic  Pressure: 
Max  Acceleration:  (g) 

Time  of  Max  Acceleration: 

Max  Heating  Rate:  (W/cmA2) 
Time  of  Max  Heating  Rate: 

Heat  Load:  (J/cmA2) 

Final  Orbit  Parameters: 
Semi-major  Axis  (km) : 
Periapsis  (km) : 

Apoapsis  (km) : 

Apoapsis  Altitude  (km)  : 
Eccentricity : 

Capture  Time: 

Final  Inclination: 

Total  Aerocapture  Mass  (kg) : 
Propellant  Mass  to  Circulariz 
Estimated  Front-shield  mass  ( 


66.52 

149.76  sec  (2.50  min) 

1354.90 

114.93  sec  (1.92  min) 

2.29 

114.93  sec  (1.92  min) 

27.67 

91.29  sec  (1.52  min) 

4612.22 

3549.02 

3408.13 

3689.92 

300.00 

0.04 

178.05  sec  (2.97  min) 

70.00  deg 

34.97 

After  Aerocapture  (kg) :  11.91 

r)  :  23.06 


Pure  Propulsive  insertion: 

Isp  =  330.000000 

Propellant  Mass  (kg):  386.401888 
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APPENDIX  F:  MINIMUM  INITIAL  FLIGHT  PATH  ANGLE  FOR 

AEROCAPTURE  AT  MARS 

Data  Summary: 

Inertial  Velocity  Components: 

Arrival  V  infinity  (m/s):  5200.00 

Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Atmospheric  Exit: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) 

Delta-V  (m/s) : 

Rotating  Velocity  Components 
Atmospheric  Entry: 

Speed  (m/ s) : 

Heading  (deg) : 

Flight  Path  Angle  (deg) : 


Atmospheric  Exit: 

Speed  (m/ s) : 

3358.39 

Heading  (deg) : 

-73.83 

Flight  Path  Angle  (deg) : 

4 .71 

Delta-V  (m/ s)  : 

3757.48 

Trajectory  Analysis: 

7115.88 

-75.27 

-10.93 


3435.29 

-69.85 

4.60 

3743.54 


7178 . 83 

-73.40 

-10.83 


Total  Pass  Time: 


352.25  sec  (5.87  min) 
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Minimum  Altitude  (km) : 

Time  to  Min  Alt: 

Max  Dynamic  Pressure  (Pa) : 

Time  of  Max  Dynamic  Pressure: 
Max  Acceleration:  (g) 

Time  of  Max  Acceleration: 

Max  Heating  Rate:  (W/cmA2) 

Time  of  Max  Heating  Rate: 

Heat  Load:  (J/cmA2) 

Final  Orbit  Parameters: 
Semi-major  Axis  (km) : 

Periapsis  (km) : 

Apoapsis  (km) : 

Apoapsis  Altitude  (km) : 
Eccentricity: 

Capture  Time: 

Final  Inclination: 

Total  Aerocapture  Mass  (kg) : 
Propellant  Mass  to  Circularize 
Estimated  Front-shield  mass  (k 


56.32 

93.33  sec  (1.56  min) 

3298.44 

78.99  sec  (1.32  min) 

5.57 

78.99  sec  (1.32  min) 

39.76 

65.53  sec  (1.09  min) 

3138.83 

3391.27 
3092.61 
3689.92 
300.00 
0.09 

98.27  sec  (1.64  min) 

70.00  deg 

42.1 

After  Aerocapture  (kg):  26.41 

) :  15.69 


Pure  Propulsive  insertion: 

Isp  =  330.000000 

Propellant  Mass  (kg):  386.401889 
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